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ABSTRACT 


The primary goal of this thesis is to identify the "state-of-the-art" in orbit-decay- 
induced uncontrolled reentry/impact prediction methods, with an emphasis on the physics 
of the final few revolutions to impact. This was accomplished through a comprehensive 
literature survey from the 1950’s to the present of unclassified military and civil 
databases. The results of the survey show that the current U.S. and international 
reentry/impact prediction methodologies are based on analysis which is over 30 years 
old. Of the various "extensions" to the current reentry theory, of which the NORAD 
method is recognized as the international standard, there does not appear to be any one 
method which is singularly superior to the others. It has also been shown that numerous 
reentry investigations made simplifying assumptions due to insufficient data needed to 
accurately model reentry and also because of computing limitations of their day. Also, 
current deterministic dynamic models appear to inadequately describe the actual 
uncontrolled reentry process, due to a lack of observational data, uncertainty in 
determining aerodynamic coefficients, atmospheric density, and point mass modeling 
where changes in vehicle configuration, attitude and lift are neglected. Stochastic and 
Statistical methods could be applied to the current methodology, to better analyze the 
various uncertainties, which could help to improve the overall predicted impact time and 
location; however, further research into these methods along with the physics of 


uncontrolled reentry is necessary. 
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I. INTRODUCTION 
A. HISTORY 

Satellite or spacecraft reentry related work has been done in this country for nearly 
five decades. With the advent of the rocket in World War II, the military saw the first 
application of spacecraft reentry prediction as the ballistic trajectory calculation of the 
target aimpoint. Prior to this, several nations had "rocket societies,” however, none had 
actually applied the use of "space" to a technology or an industry. Thus, the first 
application of space was by the military for the purpose of war. [Ref. 1:p. 13] 

On October 4, 1957, the former Soviet Union changed the perception of space held 
by most Americans and possibly the world. It was the launch of Sputnik I that caused, 
then Senator, Lyndon B. Johnson to remark: 

That sky had always been so friendly, and had brought us beautiful stars and 
moonlight and comfort; all at once it seemed to have some question marks all over 
It. [Ref. 2:p. 13] 

When the former Soviets again launched another rocket barely one month later and 
sent a dog into space, onboard Sputnik II, the world now saw the first space traveler and 
the dream of humans in space became more real than dream [Ref. 1:p. 20]. 

After a dismal failure of the United States’ first attempt to launch a satellite into 
orbit with Vanguard I, which was dubbed by some reporters as "Kaputnik” after it 


exploded on the launch pad on December 16, 1957, Explorer I was successfully launched 


on January 31, 1958. The space race had begun, and both countries were pushing to 
have the first human in space. [Ref. 1:p. 20], [Ref. 2:p. 16] 

With the goal of putting humans into space came the necessary requirement of 
providing for the safe return of those humans back to Earth, unlike the first dog in space 
which perished after about a week when the oxygen supply was exhausted [Ref. 2:p. 16]. 
This was the beginning of the most intense research period into the reentry process Over 
the entire history of space exploration. 

By April 1972, there were at least 44 reported instances where man-made space 
objects had impacted on the Earth [Ref. 3:p. 383]. By March 1978, the total count of 
man-made objects placed into Earth orbit was 10,791 [Ref. 4:p. 107]. By Aug 1991, the 
cumulative count of objects ever placed into space was 21,231 with 14,417 
decays/reentries, leaving 6,814 objects in Earth orbit [Ref. 5:p. v]. 

The return or reentry of manned spaceflights has been covered by the various 
media sources with varying degrees of intensity based upon the "newsworthiness" of the 
event and sensed public interest: 

Solar Max satellite plunges to Earth. [Ref. 6:p. 21(N)] 


Spacecraft’s study of sun ends today: Solar Max heads for fiery reentry. [Ref. 
7:p. A3] 


If you see a shooting star Dec. 8, make a fast wish for a deep cave. [Ref. 8:p. B1] 
1. Case Studies 
It is useful for the reader to understand the motivation for studying the 


reentry process. Obviously, the routine recovery of numerous U.S. and Soviet spacecraft 


must imply that the reentry process is adequately understood. In the case of controlled 
reentry and space vehicles designed for reentry, this may well be the case; however, in 


the uncontrolled reentry case, this 1s not as easy to conclude. 


a. Sputnik IV 

September 5, 1962, the first recorded impact of a man-made space object 
in the United States was satellite number, 1960 e1, Sputnik IV. This satellite reentered 
the Earth’s atmosphere over North America and many fragments of the reentering debris 
ended their orbit over the state of Wisconsin. The largest fragment of which was found 
to weigh approximately 21 pounds (9.49 kg) impacted in the city of Manitowoc, 
Wisconsin at approximately 0530 local time. [Ref. 9:p. 1] 

Sputnik IV, which had been launched on May 14, 1960, was designed 
to test life-support systems for Soviet manned space flight. On May 19, a planned 
deorbit maneuver failed and the pressure vessel separated from the cabin, leaving the 
cabin in a "lopsided" orbit. [Ref. 9:p. 2] There were a total of nine separate orbiting 
objects associated with satellite 1960 €1. The first of these objects to reenter was the last 
stage rocket body, which reentered on July 17, 1960. By July 1, 1961, six of the nine 
pieces had reentered. The payload reentered on September 5, 1962 leaving only one 
other object in orbit. [Ref. 9:p. 3] 

The U.S. Space Detection and Tracking System (SPADATS) predicted 
Sputnik IV would reenter on or about September 6, 1962. Moonwatch observers, teams 
of volunteers who attempted to observe reentry events with either the naked eye or 


telescopes, were notified of the SPADATS prediction on August 29, 1962. [Ref. 9:p. 10] 


It was based on these observations that the fragments of Sputnik IV were eventually 
located, with the exception of the largest fragment, which was found by a routine police 
patrol of the city of Manitowoc. Figure 1 shows the satellite ground trace over the last 
one-half revolution [Ref. 9:p. 20]. Figure 2 shows the ground trace of the last one 
hundred nautical miles [Ref. 9:p. 12]. Figure 3 shows the impact location of objects 
recovered in Manitowoc [Ref. 9:p. 17]. The objects found on the roof of the church 
annex, noted in Figure 3, consisted of 15 small spherules approximately 1/8 inch in 
diameter. These were about 325 feet further downrange from the initial mass which 
impacted near the intersection of Park St. and 8th St. 

There was no loss of life and the only reported property damage was the 


impact impression left in the street. 


b. Skylab 

In May 1973, the last U.S. Saturn V rocket launched Skylab into orbit 
approximately 270 nautical miles (nm) above the Earth. Skylab was the first U.S. space 
station and the centerpiece of the U.S. space program since the last U.S. moon mission. 
Skylab would remain in orbit until July 11, 1979. [Ref. 10:p. 1} 

Knowing the inevitability of decay and reentry of low Earth orbits, 
NASA contracted the Lockheed Missile and Space Company (LMSC) to investigate the 
predicted reentry and breakup of Skylab three years prior to its launch. This initial 
report predicted that Skylab would begin to breakup at an altitude of 65 nm (120 km) and 


that debris would fall 3600 nm downrange from the initial breakup. [Ref. 10:p. 2] 
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Figure 1: Sputnik IV Final Revolution 


[Ref. 9] 
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Figure 3: Sputnik IV Impact Area Ground Trace 
[Ref. 9] 


Based on the post-flight reconstruction of the reentry of Skylab, Marshall 
Space Flight Center (MFSC) was able to conclude that telemetry signals were still being 
sent from Skylab as it passed over Bermuda and Ascension Islands. Also, each of these 
tracking stations reported only one radar contact. It was therefore concluded that Skylab 
was intact and breakup had not begun although the altitude had decreased to 57 nm. [Ref. 
10:p. 5] Survivability underestimation has been typical of past analysis as will be 
discussed in Chapter III. 

Figure 4 shows the ground trace of Skylab over the last one-quarter 
revolution and denotes the impact footprint (debris dispersion area) [Ref. 10:p. 18]. 
Table 1 is a partial listing of recovered debris from Skylab and the recovery location 
[Ref. 10:p. 10]. 

Again, there was no loss of life, no human injury and no property 


damage as a result of this uncontrolled reentry. 


c. Cosmos 954 
January 24, 1978 marked the first uncontrolled reentry and Earth impact 
of a nuclear powered artificial satellite [Ref. 3:p. 384]. This was also the first case of 
"significant" property damage due to an artificial Earth satellite impact. It was reported 
that Canada spent over 11 million dollars and the United States spent nearly 3 million 
dollars in the location and recovery of radioactive debris [Ref. 3:p. 386]. Figure 5 
shows the reentry ground trace and impact dispersion area of Cosmos 954 [Ref. 11:p. 


303]. 
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Figure 4: Cosmos 954 Impact Area / Debris Dispersion 


[Ref. 10] 


Table I: 
[Ref. 10] 


SKYLAB RECOVERED DEBRIS (PARTIAL LISTING) 


a ioc 


Charred Fragments OwS 33.9S, 121.9E (In Esperance) 
Burned Material ows 33.9S, 121.9F (In Esperance) 
Aluminum 356 Casting Ows 33.7S, 122.1E (20 mi NE of 
Esperance) 
Foam Fiberglass ows 33.9S, 122.0E (9 mi E of 
Beam Section Esperance) 
H50 Tank ows 33.8S, 122.0E (9 mi NE of 
Aft End Esperance) 
1140 Tank OwS 33.9S, 122.1E (10 mi E. of 
Esperance) 
10° Steel Strip ows 33.9S, 122.3E (25 mi E of 
(H,0 Tank) Esperance) 
Heat Exchanger OWS 33.9S, 122.1E (12 mi E. of 
(H50 Cooler)’ Esperance) 
Segment of OWS 33.9S, 122.1E (11 mi E of 
Fiberglass Sphere Esperance) 
Insulation ows 33.9S, 122.1E (11 mi E of 
(Bulkhead) Esperance) 
Aluminum Gear ows 33.78, 122.5E£ (40 mi NE of 
and Housing (Urine Separator) Esperance) 
NoTank AM 33.28, 122.6E (60 mi NE of 
Esperance) 
Electronics Module AM 33.5S, 122.3E (35 mi NE of 
Esperance) 
N.Sphere AM 33.5S, 122.8E (49 mi ENE of 
Esperance in 
| Neridup area) 
Pressure Tank IU 33.2S, 122.6E (60 mi NE of 
Esperance) 
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Figure 5: Skylab Final One-Quarter Revolution 


[Ref. 11] 


The Outer Space Treaty of 1967, the Liability Convention of 1972 and 
the Registration Convention of 1976 are the three principal documents, drafted by the 
United Nations Committee on the Peaceful Uses of Outer Space and its Legal 
Subcommittee, which define ownership responsibilities of space objects [Ref. 12:p. 457). 
Although the Outer Space Treaty, Article VIII, provides in explicit terms the 
requirements that the registry state (country of origin) retains jurisdiction and control 
over its satellites while in space or on a celestial body, [Ref. 3:p. 389] there may be no 
legal requirement under the Liability Convention or the Rescue and Return Agreement 
for the former Soviet Union to reimburse either the U.S. or Canada for the cleanup cost 
of Cosmos 954 [Ref. 3:p. 387]. 

It is because of the potential for the loss of human life, significant 
property damage and legal responsibility based in international law that the reentry of 
uncontrolled artificial satellites must be further investigated. In the case of Cosmos 954, 
on the morning of reentry, the Soviets predicted reentry impact near the Aleutian Islands 
(52°N,173°W) and the North American Air Defense Command (NORAD) predicted 
reentry impact near Hawaii (19°N,156°W) [Ref. 4:p. 110]. The physics of uncontrolled 
reentry and the modeling of the reentry process must be better understood in order to 


improve the accuracy of reentry predictions. 


B. PROBLEM STATEMENT 
The probability of space objects reentering the Earth’s atmosphere and surviving 


to impact is relatively small; most artificial objects burn up in the atmosphere before they 
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impact [Ref. 13:p. 1]. Falling debris nonetheless presents a potential hazard to people 
and property. The most difficult problem associated with reentry impact prediction is 
not to determine what will survive to impact in as much as where impact will occur. 
[Ref. 5:p. v] 

The significant factors related to the prediction of time to impact and location of 
the reentry are a lack of observational data over the entire orbit and a lack of "precise" 
mathematical models which accurately describe the physical processes occurring during 
reentry. The major parameters which contribute to the uncertainties of the reentry 
prediction are: [Ref. 5:p. v] 

1. Atmospheric density variations--atmospheric density is strongly influenced by 
solar and geomagnetic activity, both of which are difficult to forecast. 

2. Aerodynamic force models--aerodynamic forces are a function of attitude, lift 
and drag coefficients, gas-surface interactions and gas dynamics such as 
continuum flow or free-flow regimes. 

3. Spacecraft attitude motion--the attitude of the object and how it changes with 
time is an important factor in estimating the aerodynamic forces encountered 
during reentry. 

4. Changes in configuration--ablation and fragmentation cause changes in 
configuration (profile area and mass/area loss) which can change the 
aerodynamic forces experienced either as an increase or decrease in net force. 

The lack of regularly spaced observational data over the entire orbit can severely 
handicap the efforts to predict reentry, especially in the final phase. Since most 


observational data is from radar tracking stations and the object is in a very low orbit, 


the time in view is of short duration as well as limited by the geographic location of the 
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tracking stations. The availability of good quality tracking data is required for accurate 
reentry predictions and could compensate for inherent deficiencies in the models. [Ref. 


5:p. vi] 


C. PURPOSE 
This research was initiated by the Air Force Space Command. The purpose of this 
thesis is threefold: 
1. Conduct a comprehensive literature survey in the area of artificial satellite 
reentry specifically, uncontrolled orbit decay and reentry into the Earth’s 
atmosphere. 


2. Describe the "state-of-the-art" of reentry/impact prediction techniques. 


3. Define critical areas of research where increased emphasis is required in order 
to improve the accuracy of the reentry prediction. 


The focus of this thesis is the physical processes of the uncontrolled reentry from 
the final few revolutions to impact. It is not the intended purpose of this thesis to 


examine in detail or focus on the following aspects of the reentry prediction problem: 


1. Atmospheric density models 
a. Solar activity and influence on reentry 
b. Geomagnetic influence on reentry 


2. Long-term orbital lifetime prediction 
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It is necessary, however, to discuss these aspects of the reentry process in order to 
completely study the physical phenomena associated with reentry. These topics are the 
subject matter of Chapter II and are considered the fundamental background material 


required for a further, more detailed study of the reentry process. 


D. METHODOLOGY 

The primary research methodology of this thesis was a comprehensive literature 
survey of military and civil aerospace data bases. The literature search of Department 
of Defense (DoD) data bases was restricted to unclassified work. The literature search 


was conducted at or through the following activities: 


1. Naval Postgraduate School, Dudley Knox Library 

2. Hanscom Air Force Base Research Library 

3. University of Colorado at Boulder, Astrophysics and Engineering libraries 
4. AFSPACECOM Astrodynamics Division (CNY) technical library 


5. United States Air Force Academy Library 


The secondary research methodology of this thesis was personal interviews with 
"experts" in the study of reentry or reentry related fields. The activities which were 


contacted or visited personally include: 


1. TRW Corporation 


2. The Aerospace Corporation 
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3. Phillips Laboratories 

4. Air Force Geophysics Laboratory (AFGL) 

5. Naval Space Surveillance Center (NAVSPASUR) 
6. Air Force Space Command (AFSPACECOM) 


7. United States Space Command (USSPACECOM) 


The primary data bases and periods over which the literature survey was conducted 


include: 
1. NASA 1940 - 1993 
2. DTIC 1950 - 1993 
3. IAA 1960 - 1993 
4. STAR 1960 - 1993 
5. DIALOG 1970 - 1993 


These data bases were searched with similar strategies in an effort to pull as many 
"original" articles as possible, while concurrently verifying the search process by 
producing numerous "duplicate" articles found in other data bases. The primary search 
terms included: reentry, reentry prediction, satellite reentry, atmosphere reentry, 


spacecraft, spacecraft reentry, uncontrolled reentry, reentry impact prediction, satellite, 


16 


descent trajectory, atmosphere drag, atmospheric density models, satellite orbit decay, 
reentry heating, reentry dynamics and reentry effects. 

In order to prevent confusion of common variables when derived in multiple 
sources, the equations and vanable nomenclature presented throughout this thesis are as 
presented in the original works, with the exception of minor changes made as noted. An 
example of such a change 1s the flight path angle, -y, presented in Chapter If. Since this 
is a survey of the literature, which dates back to the 1950’s, the authors have made a 
conscious effort to preserve the flavor of the individual works and no attempt has been 
made to standardize the nomenclature. However, this is an object for further study as 


indicated in Chapter VI. 
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II. FUNDAMENTALS OF ATMOSPHERIC REENTRY 


A. TYPES OF REENTRY 

The reentry trajectories of space vehicles can be classified into two types, either 
uncontrolled or controlled. The focus of this thesis is to investigate the uncontrolled 
reentry of satellites, however, three major types of controlled reentry, ballistic, gliding 


and skip, as shown in Figure 6, will be described for completeness [Ref. 14:p. 6]. 


1. Uncontrolled Reentry 
Uncontrolled reentry may be the result of an unrecoverable satellite subsystem 
failure or the end of the satellite’s operational life. The flight path angle is usually much 
less than one degree and lift is considered negligible. Typically, very large uncertainties 
in impact point predictions are created by decay-induced uncontrolled reentries. [Ref. 


15:p. 44] 


2. Controlled Reentry 
During a controlled reentry, the vehicle’s aerodynamic and heating loads are 
maintained within acceptable limits by controlling the effects of lift and drag forces on 
the vehicle throughout the flight. This is accomplished through a carefully designed 
space vehicle, flight trajectory and possibly a precision guidance system. Controlled 
reentry spans an aerodynamic flight regime from subsonic to Mach 25 and beyond for 


hyperbolic reentry. [Ref. 16:p. 231] 
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------- Trajectory 


Portion of trajectory over 
which analysis is applicable 








Figure 6: Types Of Reentry 
[Ref. 14] 
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a. Ballistic Reentry 
A ballistic reentry is characterized by sufficiently steep reentry angles 
where the force of lift is assumed to be negligible [Ref. 14:p. 2]. The ability to control 
the reentry velocity, flight path angle, ballistic coefficient and atmospheric properties 
determines the accuracy of ballistic reentry vehicles. Impact accuracies within the 
intended target vary from 20 km for shallow angle reentries, such as Mercury type 
vehicles, to an accuracy of a few hundred meters for a steep angle reentry of an 


intercontinental ballistic missile type vehicle. [Ref. 16:pp. 237-242] 


b. Gliding Reentry 

A gliding reentry is characterized by a glide slope rather than a reentry 
trajectory. During a gliding reentry, a vehicle such as the space shuttle creates enough 
lift to maintain a long hypersonic glide at a small flight path angle [Ref. 16:p. 242]. A 
measure of the vehicle’s lift that influences the descent path and cross range capability 
of a vehicle is the lift-to-drag (L/D) ratio [Ref. 15:p. 47]. By adjusting the vehicle L/D 
ratio or bank angle, oa, (the angle between the vehicle lift vector and the plane containing 
the vehicle position and velocity vectors), the extent of the range and cross-range can be 


controlled [Ref. 15:p. 44]. 


c. Skip Reentry 
A skip reentry is characterized by a vehicle whose L/D is greater than 
zero. Sufficient lift is produced to dominate the gravitational and centrifugal forces. If 


this lift is combined with a large enough initial angle of descent, a reentry trajectory with 
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one or more skips may be produced. As the vehicle begins to reenter, it reaches a 
minimum altitude where it begins to "pull-up" due to the lift force dominating over 
gravity. Eventually, the vehicle will exit the atmosphere at a reduced velocity. If the 
exit velocity and flight path angle are correctly controlled, the vehicle will achieve a brief 
orbital phase followed by a second reentry downrange from the first. [Ref. 16:pp. 245- 


246] 


B. MAJOR FACTORS INFLUENCING UNCONTROLLED REENTRY 

It is the focus of this thesis to describe the physical processes of the final few 
revolutions of a reentry body’s orbit and the reentry phase to impact. Reentry trajectory 
techniques differ from orbit determination techniques. The reentry phase is very 
dynamic as opposed to the exoatmospheric phase of orbital motion. Specifically, the 


reentry phase is characterized by: 


1. Rapidly decreasing altitude 
2. Rapidly increasing aerodynamic heating effects 


3. Rapidly changing aerodynamic load effects 


The reentry phase can be characterized by that portion of the trajectory prior to 
breakup and after breakup. Prior to breakup, the reentry trajectory equations of motion 
include parameters such as gravity, atmospheric density, ballistic coefficient, position and 


velocity. Solutions to the reentry equations of motion can be found using analytical, 
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semi-analytical and numerical techniques. Breakup can be described as that point in the 
trajectory where heating effects and load effects cause the reentry body to lose its 
structural integrity. 
1. Modeling Reentry Equations Of Motion 
a. Basic Equations of a Rigid Body 


The three basic equations of a rigid body are 


Cy (1) 
dt 
ST 2 (2) 
dt 
GH _vy (3) 
dt 
where 
H = reentry body’s angular momentum vector relative to its center of mass 
K = = reentry body’s linear momentum vector 
M_ = total moment force vector relative to the center of mass 
F = total force vector acting on the body 
r = position vector of the body 
Vs = _velocity vector of the center of mass 


Equations (1) and (2) are the kinematic and force equations, respectively. 


When these equations are coupled, they yield one second-order, non-linear, vector 
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differential equation which describes three-dimensional motion of the reentry body’s 


center of mass, given by 


Or _ 
Ot i GT a8gay (4) 


where g is the gravitational field and a,,,, represents the atmospheric drag acceleration. 
It is atmospheric drag that causes an artificial satellite to decay and reenter. 

Equation (3) describes the motion of a body about its center of mass, 
commonly referred to as attitude motion. The reentry body’s attitude is related to angle 
of attack, a, and is an important parameter which will be discussed in the next section. 

When the system of equations (1) through (3) are coupled along with the 
attitude parameters, it yields two second-order, non-linear, vector differential equations 
or six second-order scalar differential equations. These six differential equations 
completely describe the six-dimensional motion of and about the center of mass. This 
system of equations is commonly referred to as the six-degree-of-freedom model which 


will be described further in Chapter II. 


b. Modeling Gravity 


In equation (4) the gravitational field, g, may be written as an inverse 


square relationship 


g--He (5) 
A 
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where 
+= gravitational constant (Earth = 3.9865 x 10° km*/sec’) 
Or more accurately via the geopotential model. The geopotential model divides the Earth 
into three sets of geographically divided regions described by latitude and longitude as 
shown in Figure 7 [Ref. 17:p. 233}. The significance of the geopotential model is its 
ability to describe the gravitational field more accurately than a point mass model for the 
Earth. For example, in the region of reentry below 120 km, the first-order zonal 
harmonic J,, may attain a magnitude approaching that of atmospheric drag [Ref. 18:p. 
40]. 
c. Modeling Atmospheric Drag 

As previously mentioned! the force that causes an artificial satellite’s 
orbit to decay is atmospheric drag. The atmospheric drag acceleration vector, in 
equation (4), acts in the direction opposite of the satellite’s relative velocity vector and 


is given by [Ref. 17:p. 258] 


Vareg = po Vay (6) 
where 
Cy = Satellite drag coefficient 
A = aerodynamic effective cross-section area 
m = Satellite mass 
p = local neutral density of the atmosphere 
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Figure 7: Earth Geopotential Model 
[Ref. 17] 
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V, =I1V-w, x rl (satellite’s airspeed) 


Vv = unit vector in direction of (V - w, X r) 
w, = 2a rad/day 
V_—_s= Ssatellite’s inertial velocity vector 


The ballistic coefficient, from equation (6), 1s defined as 


Bee te (7) 


m 


Atmospheric drag dissipates the satellite’s energy which in turn causes 
a decrease in the semi-major axis. Drag affects high eccentricity orbits by gradually 
decreasing the apogee altitude, while maintaining a nearly constant perigee altitude, 
resulting in circularizing the orbit as shown in Figure 8 [Ref. 17:p. 258]. This 
contraction will continue until the satellite begins the reentry phase. Under certain 
simplifying assumptions, such as ignoring the rotation of the atmosphere, the analytical 
results describing the change per one revolution in semi-major axis and eccentricity for 


orbital decay are given by [Ref. 19:p. 230] 


3 
C 5 
Aa = ae 9 Rag sO: (8) 
m 0 i 
(1-ecosE)? 
C I 
he = aa ia l+ecosE *(cos E+e) dE (9) 
m <0 1-ecosE 
where 
a = semi-major axis 
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Figure 8: Drag Effect On High Eccentricity Orbits 
[Ref. 17] 
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(1) Modeling Atmospheric Density. The Earth’s atmosphere is 
primarily composed of nitrogen and oxygen. Solar radiation affects the dynamic 
properties of this medium by constantly changing the temperature, pressure, chemical 
constituents, particulate presence and electrical properties [Ref. 20:p. 5]. The inability 
to model the atmospheric neutral density is an integral part of the satellite reentry 
problem. Neutral density is defined as the density of the neutrally charged constituents 
of the atmosphere. The atmospheric density is not precisely known along the satellite 
path because it varies with geographic location, solar and geomagnetic conditions, 
altitude and time [Ref. 17:p. 257]. 

Atmospheric density models are categorized as either "theoretical," 
or "empirical" models. Empirical models describe the phenomena of the Earth’s 
atmosphere, based on a summary of observed data and are constructed independently of 
the laws of physics. Theoretical models apply the laws of physics. Empirical and 
theoretical models may have overlapping domains. For instance, data used to construct 
a particular empirical model may be smoothed and/or extrapolated based on theoretical 
considerations. Likewise, theoretical models must be compared to empirical data in 
order to define empirical parameters such as average sea level values of pressure and 
temperature and to establish assumption limits. [Ref. 21:pp. B-1,B-2] For example, the 
well known Jacchia models and the global dynamic models are semi-theoretical and semi- 


empirical in nature. The global dynamic models describe the physical and chemical 
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processes of the Earth’s coupled thermosphere-ionosphere system. Solar radiation and 
auroral input measurements are used in the physical model to predict the time-dependent 


density response. [Ref. 22:p. 1] 


A simplified analytical atmospheric model demonstrates the 


relationships between temperature, altitude and density 


p = pel th 7 
where 
Po = Sea level air density 
Z = altitude 
H = atmospheric scale height (RT, / g) 
To |= atmospheric temperature at sea level 
R = gas constant (air) 


Equation (10) is a simple analytic relationship where the atmospheric density decreases 
exponentially with altitude when gravity, temperature and chemical composition are 
assumed to be constant at all altitudes. This model represents a rough approximation of 
the atmospheric density. Additional physical properties and levels of sophistication can 
be incorporated into this simplified model in order to better describe the actual 
atmosphere. 

Empirical atmospheric density models use a variety of functions to 


describe the atmosphere. Some of these functions are common to most models. The 
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following list provides a brief description of these common functions: [Ref. 21:pp. B-13-- 


17) 


1. Altitude--density decreases with increasing altitude. 


2. Early afternoon bulge--maximum daytime temperatures cause the density to 
increase at satellite altitudes. 


3. Average 10.7 cm (F,9-7) flux--the average solar power per unit area at a 
frequency of 2800 MHz (A=10.7 cm) measured at various Earth locations. The 
10.7 cm flux is closely correlated to the extreme ultra-violet (EUV) radiation 
which heats up the upper atmosphere and is used as an indicator of the EUV 
flux, since EUV flux is absorbed at higher altitudes and is converted to heat. 


4. Daily 10.7 cm flux--accounts for the rapidly changing values of the Fio7 
measurement. 


5. Geomagnetic index (A,)--related to the activity of charged particles. Usually 
modeled as a correction to the atmospheric temperature. 


6. Winds--speeds up to 300 m/s have been observed in the upper atmosphere. 


Wind can significantly change the drag experienced by the satellite since drag 
iS proportional to the square of the velocity with respect to the surrounding air. 


Numerous empirical atmospheric density models have been 
developed since the launch of Sputnik I [Refs. 23-31]. Early models such as the Jacchia 
70 and 71 were derived from the analysis of satellite drag. These models identified the 
upper atmosphere as dependent on solar flux, geomagnetic index, diurnal, monthly and 
seasonal variations. A later model, the Jacchia 77 model incorporated composition data 
of nitrogen (N,) and mono-atomic oxygen (O) which was observed from satellite mass 


spectrometers. The mass spectrometer and incoherent backscatter (MSIS) model was 
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constructed using composition data derived from satellite mass spectrometers, 
accelerometers and ground based incoherent backscatter measurements. 

Several comparison studies have been conducted with various 
empirical atmospheric density models to determine their overall accuracy and efficiency. 
These studies indicate the following limitations and deficiencies of empirical models 
[Refs. 21, 32-34] 

1. Accuracies of models have remained relatively unchanged for the past two 
decades. 

2. Statistical analysis of measured satellite accelerometer density data as compared 
with atmospheric density model mean values and standard deviations, indicate 
mean value accuracies of approximately + 10% with standard deviations of 
approximately + 15%. 


3. Some models are significantly less efficient in terms of computational time. 


4. The F,;)7cm measurement does not adequately represent the complex interaction 
between the EUV flux and the thermosphere. 


5. The geomagnetic index, A,, or 3 hour K, does not necessarily represent the 
physical mechanism that causes the variation in atmospheric density. 


6. New model parameters such as the precipitation index (used as an indicator of 


magnetospheric activity), are needed to model the real physical variations in the 
atmosphere. 


Additionally, our ability to predict the solar flux and the 





geomagnetic field is usually difficult and unreliable at best. As a consequence, accurate 


forecasting of atmospheric density into the future is limited. This affects the 
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determination of the predicted orbital decay rates because of the dynamic dependence on 
the variations of the F,)7 and A, indices. [Ref. 35:p. 1] 

Considerable progress has been made over the past decade in the 
development of a dynamic atmospheric model of the coupled thermosphere-ionosphere 
system. The National Center for Atmospheric Research (NCAR) has developed a model 
as the thermosphere/ionosphere general circulation model (TIGCM). The pressure 
coordinate primitive equations of the lower atmospheric meteorology are the foundation 
of the TIGCM. TIGCM uses 20 minutes of CRAY-Y-MP 8/64 computer time per 
simulated model day to compute the prognostic thermodynamic, eastward and northward 
momentum equations and diagnostic equations of state and continuity. The model utilizes 
a 5° lat-long grid with 24 constant pressure surfaces distributed from an altitude of 95 
to 500 km. Density data collected from a 200 km altitude satellite was compared with 
the TIGCM calculated density and the results showed that the TIGCM calculations were 
within + 9-12% on a point by point basis along the satellite’s track. A new version, 
TIE-GCM, added an interactive dynamo model to calculate electrodynamic interactions 
between the thermosphere and the ionosphere. Results indicate this version model is able 
to provide improved determination of thermospheric density, especially during disturbed 
geomagnetic conditions. [Ref. 22:pp. 1-6] An operational version of this model, using 
a vector spherical harmonic (VSH) technique is under development at the University of 
Michigan. With further development, this dynamic atmospheric model may be able to 


forecast global atmospheric density values with errors less than 10%. [Ref. 36] 
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(2) Modeling Aerodynamic Coefficients. The aerodynamic coefficients, 


depend upon the following factors: 


1. The shape and dimensions of the vehicle 
2. The vehicle orientation to the on-coming air flow (angle of attack) 
3. The temperature and composition of the neutral atmosphere 


4. The gas-surface interaction phenomenon 


The flow phenomena encountered around complex shapes further complicates and varies 
this parameter in the aerodynamic regimes: free molecular flow, transitional flow, and 
continuum flow [Ref. 37:p. 5]. A Cp of 2.2 is typically used as a constant value in the 
free molecular flow region above 120 km. However, Cpis a function of angle of attack, 
shape and flow regime. Therefore, knowing the vehicle’s motion about its center of mass 
iS a critical factor in determining the aerodynamic forces acting on the vehicle. These 
forces dictate the orbit decay rate, reentry and impact location. Cp can vary from 2 in 
the free molecular flow regime to a value much less than 1 in the continuum flow 
regime. 

(3) Rarified Gas Dynamics. The pressure distribution and subsequent 
forces imparted by the near flowfield on a reentry body determines the forces and 
moments acting on the body via the aerodynamic coefficients [Ref. 20:p. 203]. There 
are basically five flow regimes which have distinguishable characteristics. These five 
flow regimes may be given quantitative definition using the Knudsen number (Kn). The 


Knudsen number is defined as the ratio of molecular mean free path to a characteristic 
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vehicle dimension, usually nose radius. The Knudsen number is indexed according to 
its reference mean free path, either after collision or in the oncoming flow. 

The boundary layer may be characterized by the Reynolds number 
(Re). Flow is said to be laminar when the viscous forces are sufficiently large to damp 
out oscillations caused by the dynamic forces. A low Reynolds number is characteristic 
of laminar flow. Turbulent flow is said to occur when the dynamic forces overcome the 
viscous forces, there is random mixing of particles and large momentum exchanges 
between fluid particles. When flow over a solid body reaches a critical Reynolds number 
the initial laminar flow transitions to turbulent flow. Turbulent flow is a critical factor 
in reentry since there is much more energy near the vehicle’s surface than in laminar 
flow conditions. In turn, under turbulent flow conditions, more heat 1s transferred to the 
surface of the vehicle. The five flow regimes are shown graphically in Figure 9 and 
described as follows: [Ref. 20:pp. 203-206] 

1. Free Molecule Regime--this region is where the molecular mean free path (A) 
is relatively large compared to a characteristic vehicle dimension such as nose 
radius. When molecules collide with a boundary layer they attain the state of 
that boundary after a single collision. This flow regime is characteristic of the 
uppermost portion of the atmosphere. 

2. Near Free Molecular Flow--frequently referred to as the "slip region" because 
gas molecules will acquire the momentum of the moving boundary only after 
several collisions. If, on the average, a molecule fails to acquire the momentum 
of the moving boundary after a single collision, then it is said to lack 


"accommodation." This lack of accommodation means that the temperature is 
a nearly discontinuous function of distance away from the moving surface. 
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Figure 9: Flow Regimes 
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3. Transition Regime--little is known about aerodynamic quantities such as lift, 
drag or heat transfer. This region is not treated very well analytically. 


4. Viscous Merged Layer Regime--essentially that region consisting of dynamically 
coupled shock-wave, boundary-layer interactions. The presence of a boundary 
layer on the wall alters the boundary conditions for the shock wave, 
simultaneously, the large pressure gradients due to the shock wave strongly alter 
the boundary-layer flow. Neither the shock wave nor the boundary layer may 
be treated as discontinuities. This is the region from approximately 110 km to 
75 km; this is where the "initial pitch over" occurs and reentry "starts." 

5. Continuous Regime--the region where classical fluid mechanics of high Reynolds 
number applies, here the shock wave and boundary layers are again treated as 
discontinuities. Often this regime is subdivided into four categories (subsonic, 


transonic, supersonic, hypersonic) with lines of demarcation established by 
Mach numbers. 


(4) Vehicle Profile Area. The profile area, or aerodynamically 
effective cross-section of the vehicle, is the area presented to the oncoming flow of the 
atmosphere and is a function of the vehicle’s attitude and configuration. A spherical 
satellite maintains a constant profile area. More complex vehicles that have various 
design shapes and deployed solar panels can have highly variable areas. The ability to 
model the area depends on the known dimensions and orientation of the vehicle. Since 
an uncontrolled satellite reentry is characterized by a critical system failure or the end 
of its operational lifetime, communications with the vehicle may have been lost. Under 
this condition, information on the satellite’s attitude will not be directly available [Ref. 
37:p. 4]. Additionally, the vehicle will experience increasing aerodynamic and thermal 
loads as the altitude decreases. These forces will act to change the vehicle configuration 


by deformation and removal of structures (mass) [Ref. 5:p. iv]. 
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(5) Vehicle Mass. Knowing the mass of a satellite accurately assumes 
a comprehensive knowledge of the vehicle. This information may be available for some 
U.S. satellites, however, for foreign systems this may pose a problem. As mentioned 
in the previous paragraph, aerodynamic and thermal loads can change the satellite mass. 
The ability to model changes in mass also assumes an accurate knowledge of the 
environment and how it affects the forces and subsequent motion of the vehicle about its 
center of mass. In order to develop a model capable of predicting this motion, the 


vehicle’s moments of inertia must be known. [Ref. 37:p. 4] 


(6) Vehicle Velocity. The relative velocity of the vehicle with respect 
to the Earth’s rotating atmosphere in equation (6) is 
V, = |V-w,xr|v (11) 
The maximum deceleration and heating rate experienced by a reentry body is a function 
of velocity. As mentioned before, wind can change the vehicle’s relative airspeed which 
can affect the drag or lift experienced by the vehicle. Determination of the vehicle’s 


velocity using doppler range rate observation information is relatively accurate. 
Uncertainties in the radial velocity may be as low as .166 m/s [Ref. 38:p. 138]. 
d. Modeling Aerodynamic Lift 
Lift causes the vehicle trajectory to follow a glide slope, skip path, or 
perturbed ballistic flight path. During hypersonic flight conditions, lift is generated by 


pressure forces on the lower surfaces at angles of attack caused by motion about the 
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center of mass. A perturbation from the nominal or zero lift flight path is the net result. 


The force of aerodynamic lift is defined as [Ref. 39:p. 25] 
L = 5P C,AV? (12) 


where 
C, = coefficient of lift 

Reentry trajectories with lift reduce the thermal and structural loads on the vehicle. 
During the final revolutions of orbit decay and reentry, lift is assumed 

or modeled to be very small or zero [Ref. 14:p. 2], [Ref. 39:p. 33]. However, an 

uncontrolled reentering satellite may generate a significant lift vector depending on the 

attitude, shape and motion of the vehicle about its center of mass. When this occurs, the 

lift will not be distributed equally about the nominal flight path. This may result in a 


deviation from the projected impact point. [Ref. 15:p. 47] 


e. Reentry Equations of Motion 
The reentry equations of motion can be derived by the mathematical 
transformation of the second-order, nonlinear, vector differential equation (4). Several 
analytical and semi-analytical theories describing a satellite’s shallow reentry equations 
of motion have been developed over the last thirty five years [Refs. 14,40,41,42]. 
With the advent of manned space exploration and recoverable probes, it became 
necessary to develop accurate theories of the entry phase, during which the 


altitude, velocity, deceleration and the heating rate vary rapidly. [Ref. 41:p. 2] 


In the derivations of these theories, several strong physical assumptions were made: 
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1. Spherical Symmetry--the Earth and its atmosphere are spherically symmetric. 


2. Non-rotating Atmosphere--the rotation of the Earth’s atmosphere which is 
approximately equal to the angular velocity of the planet is neglected. 


3. Exponential Atmosphere--the atmospheric density decreases exponentially with 
altitude. 


4. Gravitational Field--the gravitational field is assumed to be constant during 
reentry at all altitudes. 


5. Coordinate System--a non-rotating two-dimensional inertial coordinate system 
with the origin at the center of the Earth. 


(1) Fundamental Equations of Entry Dynamics. The exact reentry 
equations of motion derived by Loh, reference [40], are presented to show the basic 
relationship of the vehicle’s force vectors along the radial and normal direction to the 
flight path in an inertial coordinate system as shown in Figure 10 [Ref. 40:p. 19]. This 
presentation serves as the foundation for the development of more sophisticated theories 
to be presented in Chapter III. Physical assumptions 1,3,4, and 5 were used in the 
derivation of these equations. 

The magnitude of the relative aerodynamic velocity, V,, does not 
equal the magnitude of the inertial velocity, V, due to the reentry body’s trajectory 
through the moving atmosphere. These velocities and their geometric relationships are 


shown in Figure 11 {[Ref. 40:p. 18] and are related by the following equation 


(13) 
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Figure 10: Inertial Coordinate System 
[Ref. 40] 


Local vertical 


wae eee Ee ee 





South 


Figure 11: Relationship Between Relative And Inertial Velocity 
[Ref. 40] 
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where 


V. = 1519 cos 6 (ft/sec) (velocity of the Earth’s surface at a specific latitude) 

6 = latitude angle 

v = @ (flight path angle to the local horizon from Figure 10) 

The force components from Figure 10 are defined as 
F, = -Dcos $+Lsin ¢-mgsin p = m-< (Veos >) (14) 
F, = Leos $+Dsin $-mgcos p = m-<(Vsin ) (15) 
where 
¢ =yty 
V = V, (relative velocity) 


By resolving equations (14) and (15) in the velocity direction, replacing L and D with 
equations (12) and (6) respectively, and by further rearrangement of the terms and 


simplification, according to Loh, the exact equations of motion are 


dcosy , (1 soo (8-3). 2 (2)(Ca 16) 
dp BR) ep \V 2 \D}\ mp 
(eae We 

gh —- gk) _(2)1 a7) 
dp mB} sin y BR} p 
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where 


Ry) = radius of the Earth 
8 = = inverse atmospheric scale height (Mg/RT) 
M = mean molecular weight 


The exact equations of motion of entry dynamics cannot be solved 
analytically, however, solutions can be found using numerical methods. First-order 
approximate analytical solutions have been developed by restricting the equations to 
limited regions of application. [Ref. 43:p. 25] The ability to accurately model and solve 
these equations depends on the knowledge of the initial conditions: the initial position, 
and velocity, the vehicle’s area and mass, the neutral atmospheric density, and the 


aerodynamic coefficients. 
2. Modeling Breakup 
a. Reentry Body Structural Mechanics 
A Satellite will experience structural and thermal loads during reentry into 


the Earth’s atmosphere. The structural response of the body during reentry may depend 


on the following factors: [Ref. 44] 


1. Static load effects 
2. Dynamic load effects 


3. Thermal load effects 
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The coupled effect of these factors may determine when the body will experience 


structural failure. 


b. Modeling Reentry Heating Effects 

The reentry process 1s essentially that region of flight where the vehicle’s 
velocity into the atmosphere is reduced and its kinetic energy is converted into thermal 
energy in the surrounding medium. Since breakup is determined, in large part, by when 
the outer surface reaches its melting point, reentry heating directly affects survivability. 
This will be discussed in detail in Chapter III. The conversion fraction of kinetic energy 
to heat energy is a function of the satellite’s shape, velocity, and altitude. At very high 
altitudes (free-molecule flow region) the heat energy is developed almost entirely at the 
surface of the vehicle and up to one-half of the lost kinetic energy may be converted into 
heat in the vehicle body. [Ref. 43:p. 191] At low altitudes (continuum flow region) the 
heat energy will appear in the area between the shock wave and the body. This heat is 
transferred from the hot gas to the vehicle by conduction and convection through the 
viscous boundary layer which is adjacent to the surface of the vehicle. Radiant heating 
of the vehicle from the hot gas also occurs and this contribution to the surface heating 


is dependent upon: 


1. Bluntness of the vehicle’s leading edge 


2. Excess orbital velocity of the vehicle 
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Excess orbital velocity is defined as an appreciably higher velocity than the circular 
orbital velocity for a given altitude. [Ref. 43:p. 192] 

The rate at which the vehicle heats is not exclusively dependent upon this 
energy conversion fraction, it is also dependent upon the rate of kinetic energy loss by 
the vehicle. For the case of natural orbital decay, reentry is at small flight-path angles 
and the deceleration is very slow in the upper atmosphere. The surface heating rate 1s 
relatively low despite the high conversion fraction in this instance; therefore, the 
dominating factor is the low rate of kinetic energy loss by the vehicle. Conversely, at 
steep reentry angles, where deceleration occurs rapidly, the surface heating rate is high 
although the energy conversion fraction is low. 

The total heat input into the reentry vehicle depends upon the time of heating 
as well as the heating rate. If the energy conversion rate were constant, the total heat 
input would simply be a fixed fraction of the initial energy and the type of reentry would 
not be of significance. 

Three aspects of the aerodynamic heating process are significant, namely: 


[Ref. 40:p. 181] 


1. The total heat input, Q. 


2. The time rate and maximum time rate of local stagnation region heat input per 
unit area (dH/dt) and (dH,/dt),.,,. 


3. The time rate and maximum time rate of average heat input per unit area 
(dH,,/dt) and (dH,,/dt),.,.. 
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The time rate of average heat transfer per unit area is given by equation (18) below. The 
average heat transfer is simply the total heat input divided by the time of input. Overall, 
reentry vehicle structural integrity 1s a function of the average heat input. [Ref. 40:p. 
181} The time rate of local stagnation region heat input per unit area is given by 
equation (19). Local structural integrity is a function of local stagnation region heat 
input or the generation of local "hot spots." [Ref. 40:p. 182] The total heat input is 
given by equation (20). 

(1) Reentry Heat Input. From the three previously mentioned areas of 


concern, the following generalized aerodynamic heating equations are given [Ref. 40:pp. 
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where 
H,, = average convective heat transferred per unit area, ft-lb/ft 
H, = local stagnation region heat input per unit area 
Q = convective heat transferred, ft-lb 
and 
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= equivalent skin-friction coefficient 





velocity, ft/sec 

= atmospheric density, slugs/ft? 
= Knudsen number = )/Re 

= Reynolds number 


= molecular mean free path 


coefficient of viscosity, lb-sec/ft? 


Prandtl number, subscript e indicates "entry" value 


ratio of specific heats 


nose or leading edge radius, ft 


constant = 6.8 to 15 x 10° 


o 
21/7 


(22) 





I 


= flight path angle, positive for descent 


surface area, ft? 


coefficient of friction (local conditions) 


specific heat at constant pressure 
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Cy, = specific heat (local conditions) 
The above equations make the following simplifying assumptions: [Ref. 40:p. 181] 
1. Radiative heat transfer from the surface generally does not appreciably influence 
convective heat transfer to a vehicle; therefore, it is disregarded. 
2. Effects of gaseous imperfections may be neglected. 
3. Shock-wave boundary-layer interactions may be neglected. 
4. Prandtl number is constant. 


5. Reynolds analogy is applicable. 


If the heat transferred to the reentry vehicle is expressed as a 


fraction of the total kinetic energy, then [Ref. 20:p. 138] 





KE = 5mV; (23) 
where 
m = mass of reentry vehicle 
Ve = reentry velocity 
2 
g- —2 (24) 
mV; 
where 
= fraction of heat transferred to the reentry vehicle 
Q = total heat input 
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This may be rewritten as 





0-2 |S] fem os 
where 
Cp, = coefficient of friction 
Cy, = coefficient of drag 
S = reference area (usually nose tip) 
Sw = wetted surface area 


and where the trajectory parameter is defined by [Ref. 20:p.138] 
_ 2S ee (26) 
2msiny,  2Psiny, 
when 
Ye = feentry flight path angle 
At small flight path angles, such as reentry from orbit decay, ap, will be very large, 


therefore, equation (25) reduces to 


o-1 4 (27) 
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It is now recognizable that the fraction of the initial kinetic energy, transferred to the 
reentry vehicle by convective heating, is one half of the ratio of the friction drag to total 
drag . [Ref. 20:p. 138] 

In the case of ballistic reentry at small angles of reentry, both the 
maximum heating rate and total heat load increase as the effective mass-area ratio 
(m/C,A) or ballistic coefficient increases [Ref. 40:p. 195]. 

Assuming a constant ballistic coefficient throughout the reentry process 


yields a different relationship [Ref. 40:p. 198] 





Q- Le (28) 
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Therefore, reentry at small angles of inclination reduces the maximum heating rate of the 


vehicle, however, the total heat load is increased. [Ref. 40:p. 198] 
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(2) Reentry Heating Rate. If the heating rate per unit time is defined 


[Ref. 20:p. 135] 


. 4 32 
g A (32) 


then the total heat-transfer rate may be written [Ref. 20:p. 136] 
@ = J, dds (33) 
After some simplification it is possible to write [Ref. 20:p. 137] 


@p S, V° (34) 


g- 4 


this assumes 
| Cys = CSy, (35) 


It is now possible to define dq/dt as the average rate of change of heat transfer per unit 


area 


oF p y3 (36) 
4 





q = 
It can be shown that dq/dt is a maximum when pV? is a maximum. [Ref. 20:p. 137] 


(3) Ablation. As previously discussed, the heating experienced by 
reentry at small flight path angles is different from that of reentry at large flight path 
angles. The flight duration is much longer for the first case and even with a lower 


maximum heating rate, the total heat input exceeds that of the larger reentry angle. The 
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predominant cooling mechanism of uncontrolled reentry of vehicles not designed to 

survive reentry is ablation. Ablation is the term which generally applies when: 
...there is a removal of material (and an associated removal of heat) caused by 
aerodynamic heating, and therefore embraces, melting, sublimation, melting and 
subsequent vaporization of the liquid film, burning and depolymerization. [Ref. 
15:p. 198] 

Ablative materials are measured as "effective" depending on their 
capacity to dispose of heat (latent heat) by convection in the liquid film, and by 
convection in gaseous form in the boundary layer. Previous work has shown that when 
a material experiences a large percentage of total mass loss as vaporization, it is a more 
effective ablative material. Sublimation is the process whereby all the mass loss 
undergoes vaporization, there is no liquid film, and therefore is an excellent method for 
removing large amounts of heat. For these reasons, materials which undergo sublimation 
at reasonably high temperatures are generally more efficient at removing heat and thereby 
reducing the thermal load on the reentry vehicle. [Ref. 40:p. 204] 

Radiative cooling is another mechanism by which heat is removed 
from the reentry vehicle during reentry. Radiative cooling and ablation combined are an 


effective pair in balancing the heating effects of reentry for a lifting body. For a non- 


lifting body radiative cooling is inefficient. 


C. GENERAL DESCRIPTION OF THE REENTRY PHASE 
A decaying satellite possesses a large amount of kinetic energy due to its velocity 
and potential energy because of its altitude above the Earth’s surface. As the satellite 


encounters the atmosphere, a shock wave forms ahead of the vehicle which heats up the 
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atmosphere surrounding the vehicle. This enveloping layer of incandescent atmosphere 
causes the vehicle’s temperature or thermal load to continually increase as it penetrates 
into an increasingly denser atmosphere. During this phase, the velocity continues to 
decrease as the kinetic energy is converted into heat through the atmospheric drag. If 
all the satellite’s energy were converted into heat and contained within the vehicle, then 
there would be more than enough energy to vaporize the vehicle. However, this is not 
the case. A large part of the total energy is diverted away from the vehicle by two 
processes. The first process unloads a major fraction of the heat into the atmosphere by 
a strong shock wave mechanism. The second process involves the radiation of heat away 
from the hot surfaces. [Ref. 43:pp. 1-2] 

The vehicle also experiences structural loads, which are a combination of 
aerodynamic and inertial loads [Ref. 45:p. 5]. Atmospheric drag forces usually cause 
a reduction in the vehicle’s kinetic energy. Centrifugal and lift forces cause accelerations 
normal to the direction of the motion. 

Aerodynamic lift and drag forces vary directly with the square of the vehicle’s 
velocity, V’, and with the atmospheric density, p. Deceleration of the vehicle is a 
product of two quantities, density and velocity. As the satellite penetrates further into 
the atmosphere, the density increases rapidly resulting in a corresponding decrease in 
velocity due to drag. Initially, deceleration increases as shown in Figure 12 [Ref. 43:p. 
7}. However, at some altitude the velocity begins to decrease at a faster rate than the 
increasing density, which results in a maximum deceleration. Additionally, the 


maximum heating rate, ~ pV’, occurs at an altitude somewhat higher than the maximum 
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Figure 12: Changes During Atmospheric Reentry 


[Ref. 43] 
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deceleration [Ref. 20:p. 98], [Ref. 43:pp. 6-7]. A maximum deceleration of up to 8 
g’s, as shown in Figure 13, [Ref. 14:p. 30], occurs during this phase [Ref. 39:p. 36). 
When the satellite’s skin temperature and structural load become sufficiently high, 
the vehicle will start burning and breaking up. Solar panels and other projections such 
as antennas will separate at the earliest stage, while heavier pieces will breakup at lower 
altitudes. Based on predicted and observed data, satellite breakup commences at an 
altitude between approximately 75-120 km [Ref. 10:p. 5], [Ref. 46:p. 4], [Ref. 47:p. 
39]. The resulting debris from the breakup will impact the Earth’s surface provided it 


Survives the reentry heating process. 


D. CURRENT REENTRY THEORIES 

Predicting reentry time and impact location relies upon observational data of the 
reentry vehicle. Based upon observed (measured) position and velocity, over the orbital 
path and ideally equally distributed over that path, an algorithm is used to calculate an 
elliptical orbit which best fits the observational data. If the algorithm attempts to model 
the physical reentry process, it will be referred to as a physical model by the authors. 
Another type of model to be discussed is the King-Hele or mean motion type of model 
which neglects certain physical aspects of the reentry process and focuses more on the 
observational data. The algorithm which "fits" the observational data to an elliptical 
orbit may also be used to propagate or predict future orbital locations as a function of 


time. These types of algorithms are commonly referred to as propagators. 
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Figure 13: Maximum Deceleration vs Initial Flight Path Angle 
[Ref. 14] 
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1. Physical Modeling Theory 

The current propagator in the United States, for reentry prediction, is the 
special perturbations (SP) model. This model is maintained by the Air Force Space 
Command and is the standard for reentry prediction in the United States. 

Orbital periods of less than 87.5 minutes are defined by the Space 
Surveillance Center (SSC), Cheyenne Mountain Air Force Base, Colorado, as decaying 
orbits [Ref. 48:p. 3]. The SP model uses numerical methods to incorporate zonal, 
sectoral and tesseral orbital perturbations in the calculation of decaying orbit reentry 
predictions. Gravitational perturbations are modeled by mapping the Earth into small 
grids, which allows for enhanced resolution of the geopotential. Third body gravitational 
effects, sun, moon, and planets, may also be modeled. Third body gravitational effects 
are used primarily in the propagation of highly eccentric orbits with large apogee 
altitudes. The atmosphere is modeled using the Jacchia-Nicolet model which considers 


the following: [Ref. 48:p. 9] 


1. Diurnal bulge 
2. Solar activity 
3. Geomagnetic activity 


4. Semiannual variation 


A technique known as differential corrections is used to "fit" the observations 
of the reentry vehicle into the best orbit. The differential corrections are a mathematical 


means of determining a single orbit path (ellipse) consistent with the observed data (the 
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only "known" information in the absence of active telemetry data). The need for the use 
of differential corrections arises from the fact that the model used to propagate the orbit 
into the future has inherent deficiencies. These deficiencies are then reflected as errors 
in the prediction of where the vehicle is supposed to be in its orbit as compared to where 
it is observed to be. Through a series of fitting observations to an orbit and updating the 
predicted orbital path, the reentry of the vehicle is calculated as a “time” when the 
altitude will reach a specified minimum value. This lower limit value of altitude is a 
function of the atmospheric density model. The impact dispersion area or footprint is 
calculated as + 15 minutes of that reentry time. The sub-satellite ground trace is used 
to describe the impact area on the Earth’s surface. Table 2 shows the format of 
observational data as used in the differential corrections process [Ref. 49]. Figure 14 
shows the differential corrections display and Figure 15 shows the Tracking and Impact 

Prediction (TIP) display used at Cheyenne Mountain AFB. [Ref. 49] 
It now becomes obvious that there are two extremely significant issues at 

hand: 

1. The +15 minute window equates to approximately 1/3 of one complete 

revolution of the reentry vehicle’s orbit in the decay phase. 


2. The lack of observational data in the decay phase may significantly bias the 
predicted reentry time. 


58 


OBSERVATIONAL DATA FORMAT 


Table Il 


[Ref. 49] 
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The accuracy of the reentry prediction is directly limited by the ability to 
observe the reentry process as well as indirectly by the inherent deficiencies in the 


models used to represent the physical reentry process. 


2. Mean Motion Theory 
In the early development of his mean motion theory of satellite lifetime 
prediction, King-Hele defined orbital lifetime as the time remaining until the eccentricity 
of the orbit reached zero [Refs. 19, 50-55]. This was a prediction of reentry based on 
the circularization, or contraction of the orbit, due to atmospheric drag. A principal 
assumption in this theory was that perigee height remained constant. When oscillations 
in perigee height were considered (which occurs when an oblate atmosphere is modeled), 
the theory had to be revised in order to take into account a zero eccentricity prior to 
reentry [Ref. 54]. The definition of orbital lifetime was then redefined in terms of mean 
motion, n, where a value of n (chosen by the user) represents an "end value" of a 
circular height. For example, a value of n equal to 16.5 rev/day would correspond to 

a circular orbit altitude of 150 km. 
In the prediction of orbital lifetime, the observed rate of change of orbital 
period, dT/dt, was the entering argument (where T is the orbital period) [Refs. 19, 50- 
53]. This was equivalent to using the rate of change of mean motion, dn/dt, since the 


two are related by [Ref. 54:p. 5] 


(37) 


EY 
n=— 
T 
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where 
n = rev/day 
T = days 
In its simplest form, orbital lifetime is now expressed in terms of the rate of change of 


mean motion as 


fee 2 (38) 
7 


i 
l 


lifetime in days 
dn/dt = rev/day’ 


Q 


and expanded by [Ref. 54:p. 8] 


function of e, n and H (density scale height) 


Z 3enl 


o(Z) 7, (39) 
Q = — a 11-0. 1] 





! 


where 
I,,1, = Bessel functions of the first kind, of degree 0 and 1, with argument z 
J = 0.3€, - 0.025 
Z = ae/H 


For the circular orbit, the lifetime may be written [Ref. 54:p. 9] 
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_ 3H_n(l-e™") (40) 
2an 


L 


where 
H,, = scale height H km below initial altitude 
h = decrease in altitude until reentry (120 km) 

In this form, density scale height is the only parameter which is not directly 
derived from orbital data (observations). It can be shown that H is less sensitive to solar 
activity than a neutral density atmosphere model, therefore, uncertainties in solar activity 
will have less influence on the orbital lifetime prediction [Ref. 54:pp. 9-10]. Perhaps 


more importantly, the mean motion technique does not require any knowledge of: 


1. Satellite attitude 
2. Satellite size 

3. Satellite shape 
4. Satellite mass 


5. Aerodynamic characteristics 


These factors are all accounted for in the observation of mean motion and rate 
of change of mean motion. In summary, the mean motion theory is dependent upon five 


parameters: [Ref. 54:p. 12] 


1. mean motion (n) 

2. eccentricity (e) 

3. daily solar activity index (F197) 

4. average solar activity index (Fj97) 


5. geomagnetic index (A,) 


The last three terms above are incorporated into the calculation of atmospheric scale 


height. 
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HI. FORMULATIONS AND SOLUTIONS OF REENTRY 


This chapter describes the "state-of-the-art" formulations and solutions (analytical, 
semi-analytical, and numerical) of decay-induced reentry as determined from the 
literature survey. Fundamental physical processes that were introduced in Chapter II, 
including reentry equations of motion, rarefied gas dynamics, reentry heating and reentry 


body breakup are presented. 


A. ANALYTICAL REENTRY EQUATIONS OF MOTION 


1. Chapman’s Approximate Analytical Entry Equations of Motion 
Chapman derived a nonlinear second-order differential equation by reducing 
equations (16) and (17) and by introducing a set of completely nondimensionalized 
variables [Ref. 14:pp. 3-14]. In the development of this equation, physical assumptions 
1 through 3 from Chapter IJ were used. Additionally, two mathematical assumptions 
were used: 

1. The fractional change from the center of the planet, dr/r, in a given increment 
of time is much smaller as compared to the fractional change in velocity,d(V 
cos y)/V cosy, given by the mathematical expression | d(V cos )/V cos y | 
> |dr/r|. 

2. The flight path angle, -y, related to the local horizon direction for lifting vehicles 
is sufficiently small such that the lift component in the horizontal direction is 


small as compared to the drag component, given by the mathematical 
expression, 1 > | (L/D) tany|. 
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It is erroneous to think that these assumptions will restrict Chapman’s analysis to 
entry trajectories with small flight path angles and small lift-to-drag ratios as many 
authors have believed. On the contrary, these assumptions applied simultaneously, 
constitute a well balanced set of hypotheses and make Chapman’s theory applicable 
to a large family of entry trajectories. [Ref. 42:p. 179] 


As previously mentioned, Chapman introduced a set of dimensionless 


(independent and dependent) variables. The independent variable is 


— Vcosy _ U 





= id (41) 
and the dependent variable is 
ee | TG (42) 


Zi | 6B 


Chapman’s final derivation is given by 


ieZ weal = Q =~) ogg 5 | - Vir & cos >| fe) 
u uZ D 


l. . Sh 4, 
where 8 = Mg/RT and the numbers associated with the brackets in equation (43) 


represent the following physical quantities: 


1. Vertical acceleration. 
2. Vertical component of the drag force. 
3. Force of gravity minus the centrifugal force. 


4. Lift force. 
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For the specific case of the shallow reentry of a satellite from a decaying 
orbit where the flight path angle is very small, the Z function equation (43), can be 


written in the following form 





ud (2 -2|- 02) . reo (44) 


du\du su 


with the initial conditions of 


uu=1, y,=0 (45) 


and where the corresponding boundary conditions for decaying orbits are 


Zl) =0 , Z’(1) =0 (46) 
Integrating equation (44) by numerical techniques for each designated initial 
velocity, initial flight path angle and L/D ratio generates a solution. The results of the 
solution are applicable to any vehicle of arbitrary dimension, size, or mass. Figure 16 
shows the solution graph for a non-lifting vehicle (L/D=0) [Ref. 14]. Figure 17 shows 


the solution graph of various L/D values using the parameter [Ref. 14] 
Ver |= (47) 
ae 
mu 
The ratio u/u, for the ordinate from Figures 16 and 17 is equal to equation (41). 


Additionally, the following engineering quantities of interest during the 


reentry process can be calculated from the Z function: 
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Figure 16: Z Function Solution Graph For L/D=0 
[Ref. 14] 
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Figure 17: Z Function Solution Graph For Various L/D’s 


[Ref. 14] 
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1. Deceleration 

2. Descent angle 

3. Range 

4. Time 

5. Density-velocity relationship 
6. Dynamic pressure 

7. Reynolds number 

8. Heating rate 


9. Total heat absorbed 


2. Loh’s Second Order Unified Solution of Entry Dynamics 

Loh derived a general second-order solution of reentry mechanics that covers 
the entire range of initial flight path angles and L/D ratios [Ref. 40:pp. 25-36]. Figure 
18 shows the entire range of the second-order solution as compared with several other 
analytical approximate solutions available in 1963, including Chapman’s approximate 
analytical theory [Ref. 40:p. 26]. 

The second-order unified solution can be derived from equations (16) and 
(17) by using approximations, a binomial series expansion, and integration techniques 


for a small flight path angle,y, and is given by the following equation 
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Figure 18: Loh’s Second-Order Range Of Applicability 
[Ref. 40] 
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y? Vi; 49 
2 | irr 5) Xr = =a ( ) 
BR, BR, 
Ry = radius of the Earth . 
f = condition at beginning of unpowered glide or ballistic entry 
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: (50) 
(Cy (Gy | 1,1 2, 
| 720C; 240C/ | 120C/ 12¢; oe, oe 
-] 

C R 
c= - lL] |S] -(1)_cosy_ | &-; (51) 

21D] | mB R, ie Vv? 

p pee ee 
m 


For the specific case when L/D=0, equation (48) can be simplified and 


rewritten as 


WS 


in | 18% | p [1-2 psi. (52) 
| 4 mB p COS ¥,-COS ¥ 


Figure 19 shows a comparison between Chapman’s first-order solution and 
Loh’s second-order solution for a small L/D=0.1 [Ref. 40:p. 63]. Both solutions offer 
about the same degree of accuracy. 
However, the first order solution is limited by the conditions that (L/D) tan + must 
be smaller than 1 and the initial angle of inclination -y; must be very small; the 
second order solution does not suffer these particular limitations. When L/D=0 
and +; = 0, where the first order solution was not available previously, the second 


order solution offers for the first time a satisfactory solution in this region.[Ref. 
40:pp. 50-51] 


3. Yaroshevskii’s Entry Theory 
Yaroshevskii developed a semi-analytical reentry trajectory theory which was 
originally published in a Soviet journal, Kosmicheskie Issledovaniya, in 1964 [Ref. 42:pp. 
158-176]. [An English translation of this article was not available to the authors. ] 
Using some simplifying assumptions, he derived a nonlinear second-order 
differential equation which can be integrated analytically by using series 
expansions. To some extent, Yaroshevskii’s theory is a special case of a more 
sophisticated theory developed by Chapman. [Ref. 42:p. 158] 
In the development of the theory, physical assumptions 1 through 5 apply to 
the basic reentry equations (16) and (17), from which the differential equation was 


derived. Also, for a constant angle of attack, a, the drag coefficient, Cp, and the lift 


coefficient, C,, were assumed to be a function of the Mach number. 
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Figure 19: Loh’s Second-Order vs Chapman’s First-Order 


[Ref. 40] 
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Yaroshevskii used an independent variable, x, and a dependent variable, y, 











defined as 
re -| Bovey (53) 
1 CMV 
and 
_ ©,(1)A | 7% (54) 
i a B p 
where 
To = radius of the planet 
V 
Wee (55) 
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By defining the differential relations dx and dy, using equations (53) and (54) along with 
the equations (16) and (17) and by eliminating the flight path angle, +, the second-order, 


nonlinear differential equation is given by 


1 





ae | 
dy _ _ fe GLVO)] vx) (56) 
dx? - Cod) y 


Solutions for equation (56) can be obtained by numerical integration. When 
Cy, and C, are independent of the Mach number, solutions can be obtained by integrating 


a selected series expansion of equation (54), depending on the type of reentry trajectory. 
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For the specific case of the shallow reentry from orbit decay, where 


C,/C,=0, equation (56) can be derived as 


oe : a (57) 





with the initial conditions of 
mero, y0)=0 , 2 =0 (58) 


By taking into account a singularity at y=0O, the series solution to equation (57) is 


eam [$ x7@ra2 aaa ti) (59) 
where 

a =1 

ay = 1/6 

a> = 1/24 

a, = 47/4752 


The coefficients a, can be calculated by the recurrence formula 


k-l 


2* ] 
ee 2m+1)(2m +3 : 
oe Gri 3g ONO Fi a 


‘ 1 4 (2k +1)(2k +3) 
3 


As in the case of Chapman’s theory, several engineering quantities of interest 


during the reentry process can be calculated from equation (56): 


1 


1. Time of flight 
2. Range 

3. Deceleration 
4. Heating rate 


5. Heat absorbed 


4. Universal Equations for Orbit Decay and Reentry 
Longuski and Vinh derived a set of universal entry equations of motion for 
all regimes of atmospheric flight: from the free-molecule flow regime to the near free- 
molecule flow regime where orbital motion is perturbed by air drag, through the 
transition regime to the continuum flow regime where the dynamic phase of reentry 
occurs, and to the point of impact on the planet’s surface [Ref. 41]. 
Rigorous mathematical techniques, such as Poincare’’s method of small 
parameters, and Lagrange’s expansion are applied to obtain a highly accurate, 
purely analytical theory for orbit contraction and ballistic entry into planetary 
atmospheres. [Ref. 41:p.v] 
Figures 20 and 21 [Ref. 41:pp.16-17] describe the inertial coordinate system, 


nomenclature and the aerodynamic forces of the equations of motion for a vehicle with 


a C,/Cp ratio, defined by the following equations [Ref. 41:p.10] 


= = Vsiny (61) 
ag _ Vcos ycos ¥ (62) 
at rcos @ 
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Figure 20: Coordinate System And Nomenclature 
[Ref. 41] 
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Figure 21: Aerodynamic Force Diagram 
[Ref. 41] 
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do _ Vcos ysin y 


F (63) 
t r 
AC,v? 
o = 2 gsin Y (64) 
m 
d pAC,V*cos ao y2 
yal = ——t __ -(g- cos aie 
dt 2m r 
AC, V? 2 


The exact universal equations of motion for entry trajectories for a vehicle 
inside a rotating atmosphere can be derived from a transformation of equations (61) 


through (66) using the modified Chapman variables 


_ V’cos’ y (67) 
gr 


z= — le (68) 
m 


and nondimensional independent variable 





s=[) | er (69) 


The exact universal equations for entry trajectories are [Ref. 41:pp. 11-12] 


81 


ae pr -ld_ 1,14 Ztan (70) 
ds Bo dr 2Br 2? dr 

du © _2VBr Zu | 14 cosy tany + SY (71) 
ds COS ¥Y ee 2V6r Z 














dy _ Ver Z 2 eed COS 24) (72) 


Y C 
da Se Br Z sina pel heey (73) 
ds tani cos’ D 
dQ _ vbr Zsina | % Eats (74) 
ds sini cosy 5 
di _ yBr Z cosa Cc, oe (75) 
ds cos’ c 
where 
0 = longitude of the ascending node of the osculating plane 
at = angle between the ascending node and the position vector 
1 = inclination of the orbit 
0 = bank angle 
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According to the authors, the universal equations have three advantages: [Ref. 


41:p.128] 


I= Tey are exact. 
2. They are free of any restrictive assumptions. 
3. They contain the modified Chapman Z variable, which permits a single 


trajectory solution for a specified initial velocity and flight path angle that 
applies to any vehicle of arbitrary area, mass, or Cp 


As previously mentioned, equations (70) through (75) can be used during all 
phases of aerodynamic flight and are most useful for analyzing the last few revolutions 
and the reentry phase. The accuracy of the equations depends on the readjustment of 
the value of the inverse atmospheric scale height, 6, for each layer of the atmosphere. 
[Ref. 41:p.15] 

For the specific case of reentry from a circular orbit, where C,=0, Longuski 
and Vinh derived a separate analytical theory from the exact universal equations for entry 
trajectories, due to the fact that: 

...1t does not seem possible to have a single analytic solution which is uniformly 
valid for all values of initial flight path angles because of the nature of the 
problem. In the case of atmospheric entry from circular orbit, the magnitude of 
the flight path angle, initially zero, rapidly increases, approaching 90° as the 
velocity becomes small. On the other hand, for steep angle entry, the flight path 
angle changes very little--of the order of tenths of a degree--as the nondimensional 
velocity decreases from unity to one tenth of the original value (between Mach 2 
and 3). [Ref. 41:p. 85] 


Under the condition of reentry from a circular orbit, equation (73) is equal 


to one and equations (74) and (75) are equal to zero. By using this fact and the variable 
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eee AGP a (76) 


gr cos’ 


equations (70) through (73) can be written as 


az -~BrZ tany (77) 
da 

dv _ _2ybr Zv 1- sin 1-2 (78) 
da COS y 2Ver Z Vv 

OY Stes (79) 
da v 


By dividing equations (77) and (79) by equation (78) to form a new set of equations, and 


defining the following change of variables for substitution into the new set 


Y=2Z (80) 

@ = -Ver siny (81) 

X = -Inv (82) 
1 

fae 83 

€ a (83) 


and then by expanding this new set of equations in one term, in e, and using Poincare ~’s 


method of small parameters for integration of a system with the assumed solution form 
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b= + €, (84) 


Y = ¥,+eY, (85) 


two systems of two first-order differential equations are formed [Ref. 41:p. 89] 








me. €*-! (87) 
dX Y, 

9a 
Bee @ +—° (2e%-1) (88) 
dX Y, 
ee x" | Fo (2eX-1) -42- 21 (89) 
dx ie | Y, mote 


Equations (86) and (87) can written in the following second-order differential equation 


form 
2 
oe eX 1 (90) 
dx? ve 


and with the initial conditions for the case of the shallow satellite reentry 
¥(0) = 0 , #,(0) = 0 (91) 


a series solution can be obtained in the following form 
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3 2 3 4 
2 Xana ras x AL x lll x + 20021 x (92) 
V3 314 6] 4 594 | 4 605880 | 4 


sos SE ART a EE) eel] 
Yaroshevskii’s approach was used to find the first term of the series in equations (92) 
and (93). The same approach can be applied to equations (88) and (89) to find a 
solution. 

Figures 22, 23, and 24 represent various solutions graphs for equations (92) 
and (93). The dashed line indicates the exact numerical solution while the solid line 
represents the analytical solution. Figure 22 shows the variations of aerodynamic 
deceleration, G’s (g’s = number X gravitational acceleration at a radial distance r), as 
a function of the dimensionless velocity, v, at several initial flight path angles,y,, for 
equation (92) [Ref. 41:p.120]. Figure 23 shows the variations of In ( Z / Z,) = (1- r)/H 
= drop in altitude, in units of scale height, as a function of (v) at several (y, ), for 
equation (92) [Ref. 41:p. 122]. Figure 24 shows the variation of the negative flight path 
angle, -y, aS a function of (v) at several (y; ), for equation (93) [Ref. 41:p. 121]. 

5. Attitude Dynamics of Uncontrolled Motion During Reentry 
In the previous section several analytical theories were presented describing 


the reentry equations of motion and their solutions. Strong physical assumptions were 


made in the derivations of these theories in order to describe the trajectory of the body’s 
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Figure 22: Variations Of -y vs The Nondimensional Velocity (v) 
[Ref. 41] 
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Figure 23: Variations Of In(Z/Z,) vs The Nondimensional Velocity (v) 
[Ref. 41] 
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Figure 24: Variations Of G vs The Nondimensional Velocity (v) 
[Ref. 41} 


89 


center of mass or point mass. Specifically, Chapman’s and Longuski and Vinh’s theories 
used the Z variable which permits a single trajectory solution for a specified initial 
velocity and flight path angle that applies to any vehicle of arbitrary area, mass, or Cp. 
The effect of the uncontrolled motion of a body about its center of mass on the reentry 
trajectory was not investigated in these theories. This section will examine three 
analytical investigations presented by two Russian authors in this specific area. 

Duzmak, in 1970, presented the first systematized explanation on the problem 
of the attitude dynamics of uncontrolled reentry body motion. This problem was the least 
developed in comparison to center of mass or point mass trajectories and aerodynamic 
heating [Ref. 56:p.2] 


Unifying papers on the dynamics of the motion relative to the center of mass have 
not appeared up to the present time. [Ref. 56:p. 5] 


The primary focus was the investigation and derivation of the relationship 
between a reentry body’s parameters outside the atmosphere with the body’s parameters 
in the dense layer of the atmosphere. Changes in the state of a reentry body’s motion 
relative to its center of mass during the motion along its trajectory were also investigated. 
An asymptotic approach was used on the equations of motion to solve the problem, 
specifically for the cases where the characteristic time of the entry body motion relative 
to its center of mass 1s much less than the characteristic time of motion of the center of 
mass. Additionally, a refined asymptotic method that has a significantly wider range 
of applicability was developed for those cases in which the above condition was not 


fulfilled. | This method is based on the coupling of the numerical and asymptotic 
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solutions of the equations of motion. Two-dimensional motion without restrictions on 
the shape of the body and three-dimensional motion of a body that possessed 
aerodynamic and dynamic axial symmetry were assumed in this investigation. [Ref. 
56:pp. 1-3] 

During the orbital or exoatmospheric phase, the external moments that 
produce the reentry body’s perturbed motion about its center of mass are determined by 
the laws of motion of a ngid body as described by Euler. The magnitude and direction 
of the initial angular moment vector, H, determines the state or nature of this motion. 
For instance, if initial angular moment vector during the orbital phase produces two- 
dimensional motion, then the body rotates at a constant angular velocity, w, about its Z 
or transverse axis. [Ref. 56:pp. 5-6] Figure 25 is an attitude dynamics coordinate system 
that shows the direction of the Z-axis relative to the body’s orbit [Ref. 57:p. 113]. 

Regardless of the nature of the motion, the entry body’s angle of attack, a, upon 
approaching the atmosphere can have absolutely any value at all ... Thus any 
complete solution of the problem of atmospheric entry can be obtained upon the 
necessary condition of the absence of any restrictions on the size of the angle of 
attack. One can say that the presence of large angle of attack is one of the main 
distinctive features of the problem of atmospheric entry. [Ref. 56:p. 6] 

As a body with perturbed motion relative to its center of mass enters the 
reentry or atmospheric phase, it begins to experience a stabilizing effect that is 
proportional to the increase in atmospheric density. This stabilizing effect is a property 
of dynamical systems with variable parameters where if the stiffness of the system 


increases during the perturbed motion, then the vibrations dampen out. This dampening 


effect was described, for the case of angularly misaligned missiles during reentry, by 
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Figure 25: Attitude Dynamics Coordinate System 
[Ref. 57] 
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Allen in 1957. [Ref. 58] The system stiffness effect during reentry is defined in the 


quantity 
ad (94) 
I 
where 
M,* = the derivative of the dimensional aerodynamic moment with respect to 
the angle of attack 
I, | = reentry body’s moment of inertia relative to the Z-axis 


Based on these facts, Duzmak states: 
The value of M,* increases by several orders of magnitude upon the descent 
because of the increase in density and the indicated effect of the variation in the 
system’s parameters is the main factor determining the damping of the oscillation 
... Investigating the disturbed motion upon atmospheric entry permit one conclude 


that a descending entry body represents a significantly nonlinear mechanical 
system. [Ref. 56:p. 7] 


a. Equations of Perturbed Motion 
The basic development of the equations of perturbed motion of a body 
about its center of mass during reentry is presented due to the uniqueness of Duzmak’s 
work [Ref. 56:pp. 12-28]. Equations (1) through (3) describe the system of equations 
of motion for this problem. 
The characteristic time intervals for the reentry body’s motion relative 


to its center of mass and the motion of its center of mass are, respectively, defined as 


(95) 


25 


t= Ar (96) 


characteristic rotational velocity with respect to the center of mass 


characteristic variation in the radius vector of the center of mass 


: V 
where 
0) 
Ar 
V 


characteristic velocity of its center of mass 


In the development of the equations of perturbed motion, the following 


factors and assumptions were taken into account: 


The interaction between the reentry body’s motion with respect to its center d 
mass and the motion of the center of mass is neglected. 


The hypothesis is satisfied in the upper atmosphere portion of the descent 
trajectory until the angles of attack becomes less than one radian because of the 
atmosphere’s stabilizing effect. 


The hypothesis is satisfied in the upper atmosphere portion of the descent 
trajectory since the body’s kinetic energy of the center of mass is many times 
larger than the body’s kinetic energy due to its motion relative to the center d 
mass. Due to this fact, the aerodynamic moments begin to affect the motion 
about the center of mass by stabilizing the body at significantly higher altitudes 
than the aerodynamic forces acting on the center of mass. 


The center of mass trajectory parameters are known as a function of time. 
The flight trajectory is two-dimensional. 

The average rotation of the velocity vector of the reentry body’s center of mass 
which is the curvature of the average trajectory is neglected. This effect is 
caused by the gravitational force and the body’s rotational velocity about its 
center of mass. 

The hypothesis is clearly satisfied for the sections BC and DF as shown in 


Figure 26 [Ref. 56:p.16], because the trajectory is nearly linear due to flight 
velocities of several kilometers per second. In the sections AB and FG, where 
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Figure 26: Trajectory Of The Reentry Body 
[Ref. 56] 95 


the trajectory inclination angle, 6, varies significantly, the hypothesis 1s usually 
satisfied. The hypothesis is not satisfied in the exoatmospheric section CD. 


The variation in the velocity vector orientation caused by the effect of the 
aerodynamic lift force is taken into account. 


The reentry body coordinate system, shown in Figure 27 [Ref. 56:p. 23], 


has both dynamic and aerodynamic symmetry and is used for the derivation of the three- 


dimensional perturbed equations of motion: 


10. 


The OXYZ right-handed coordinate system is fixed in both inertial space and 
relative to the reentry body. 


The body’s center of mass is denoted by, 0. 
The OXc axis is the velocity vector of the center of mass. 


The angle, a, between the OX and OXc axis 1s not the three-dimensional angle 
of attack or nutation angle. 


The plane containing the OXc, OX and OY axis is the angle of attack plane. 


The motion of the body in the angle of attack plane is defined by dx/dt and is 
directed along the OZ axis. 


The rotation of the angle of attack plane with respect to the velocity vector, 
OXc, is determined with the help of the precessional velocity, A, which is 
directed along the OXc axis. 


The body’s rotation with respect to the angle of attack plane is determined using 
the intrinsic rotational velocity, 4, directed along the OX axis. 


The angle, , between the trajectory plane and the angle of attack plane is 
defined by dy/dt = X. 


The angles a and y have the following ranges: 0 < a < 180°;0 < y < @, 


96 


Trajectory of the 
cei iadlaiahadel of mass 


\ ” 





Trace of the 

angle of 

attack's plane 
dex 
at 





E27 (E) 






XC 
Trace of the 


trajyectory's 
plane 


Figure 27: Reentry Body Coordinate System 
[Ref. 56] 
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This coordinate system can also be used to describe two-dimensional motion about the 
center of mass. 

The angular velocity vector, w, of the reentry body is the resultant of 
three rotations: da/dt; 4; and ». As previously mentioned, the aerodynamic lift force 
affects the body’s velocity vector orientation. Lift acts in the angle of attack plane, 
creating a motion of the body’s velocity vector that is in the same plane with the change 


of the angular velocity along the Z axis and is defined as 


bu, = ea (97) 
where 
L = force of lift from equation (12) 
m = mass of the body 
€ = t,/t, (small parameter) 


The reentry body’s angular velocity vector X, Y, Z components are 


defined as 
w = ptr cosa (98) 
w, = “Ah sina (99) 
wy Sg (100) 
7 dt mV 
The reentry body’s angular momentum vector, H, XYZ components are 
defined as 
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(101) 


H, = Lw, (102) 
H, = 1, (103) 
where 


I.y2 = reentry body’s moments of inertia along the X, Y, Z directions 
The main moment force, M, acting on the reentry body consists of a 
restoring moment, M,(7,a), and a small damping moment that can be projected in the X, 


Y, Z directions, is given by 











M(7,a) = m,qAL (104) 
= 2 
M,x(r,a) = m,x ac (105) 
= 2 
Myy(r,0) = myy ie (106) 
= 2 
M,2(r,a) = m,z tak (107) 
where 
ee, mex , my y, m:? z = dimensionless aerodynamic (108) 
coefficients 
q = pV?/2 
T = et 


oo 


By projecting dH/dt = M from equation (3) on the OX YZ coordinate 
system and taking equations (101) through (107) into account, the reentry body’s 


equations of perturbed motion can be written in the form 








W, 7” 
=F +1 0,0, ya, = €M, x(7,a) ow, (109) 
ee ey a 
z dt yp WW, 2, WY, eM, y 5A), 
dw, ; (111) 
I , +1 wo, -Iw,0, =M (7,0) +eM, 2(7,a) w, 
where 
w, = dAcOSa 
Wi = wy 
Ww, = WW, 


Equations (109) through (111) can be rewritten by substituting in equations (99) and 


(100) for w,, w,, and w, and by substituting w, from above, into the following form 


a tef(r,0)r = 0 (112) 
dh , (2 cosa -r) - a f(t.) - ear (113) 
dt Sin @ ae sin a 

Z 
we -\’ sina cosa+rad sina +M(7,a) tef (1,0) = 0 (114) 
dt? dt 
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where 





r= Le, (115) 
Il 
M(1,0) = a (116) 
M° 
f(7,a) = pe alba (117) 
f (7,0) = _ My W70) , L cosa (118) 
L mV sina 
S(7,a) = Bee) moe (119) 
] mV 
p2 = ob (120) 
0a 


Finally, the system of equations that describes the two-dimensional 


motion of the reentry body’s center of mass is given by 





A 
a _ |G. cord Br og g{ 1-2 2p 
dt mV V Rg, 
C.gA 
Ad = € | - =e ind (122) 
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dH 


dt 


dL 


dt 


where 


Ls a) 


Tr A © 


eV sind (123) 


i V cosé (124) 


gravitational acceleration of the Earth 


angle between the tangent to the trajectory and the local horizontal 


altitude of the flight 


R, + H 


range of flight measured with respect to the Earth’s surface 


The investigation of equations (112) through (114) in the case of known 


solutions for equations (121) through (124) is the principal content of Duzmak’s work. 


Specifically, the following major areas were investigated: 


The distinctive features of the two-dimensional motion that explains the 
interaction of the nonlinear effect with the influence of variability of parameters. 


The distinctive features of the motion that results from the transition through the 
transonic velocity range. 


The relationship of the angular momentum components with the parameters of 
motion in a increasingly denser atmosphere. This analysis clarifies the effect 
of the shape of the instantaneous characteristics, stability margin and damping. 


The effect of the reentry body’s motion parameters above the atmosphere on its 
motion in the denser layers of the atmosphere. 


The case of sinusoidal dependence of the longitudinal moment on the angle of 
attack for three-dimensional motion. 
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tt 


6. The body’s motion relative to its center of mass for a small angle of attack 
reentry (Shallow angle reentry from a decaying orbit). 


b. Effect of Motion Relative to the Center of Mass on the Motion of the 


Center of Mass 


The coupled effect of the perturbed motion about the center of mass with 
the trajectory of the center of mass occurs because of the dependence of the aerodynamic 
coefficients on the angle of attack. Duzmak neglects this interaction in the development 
of the asymptotic solutions. This assumption is based on the fact that the atmosphere 
will start to influence the motion relative to the center of mass at higher altitudes than 
the atmosphere will start to effect the motion of the center of mass. Generally, for 
perturbed motion about the center of mass in a free-molecular regime, the effect of the 
lift force is small, since it continually acts in different directions. In the denser layer of 
the atmosphere, lift can also provide some additional damping of the angle of attack 
oscillation. [Ref. 56:p. 29] 

The main effect on the trajectory is exerted by ag,,,, equation (6). For 
perturbed motion at high altitudes where atmospheric density is low and C)(qa) varies 
significantly, ag... has little effect on the trajectory. In many cases at lower altitudes, 
perturbed motion has a weak effect on the trajectory due to the fact that a reentry body 
iS: 

...able up to this instance to stabilize itself under the action of aerodynamic 


moments so that the variation in Cp in the case of perturbed motion becomes 
insignificant. [Ref. 56:p. 29] 
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For the case of hovering or similar motion where the oscillations do not 
dampen, the oscillations of the angle of attack can occur with significant amplitude. 
Under these conditions, it is essential to take into account the angle of attack oscillation 
effect on the trajectory of the center of mass. [Ref. 56:p. 317] 

The determination of the range component of the scattering landing points is 
important. This scattering is caused by the oscillations, the angle of attack and the 
variation of the coefficient of drag associated with the oscillations. This quantity 
cannot in general be determined without taking into account the interaction of 
motion relative to the center of mass and the motion of the body center of mass. 
[Ref. 56:pp. 317-318] 

In the derivation of a simple approximate method to solve this problem, 
the descent trajectory was divided into two sections: H = UH, (trajectory in the 
atmosphere’s tenuous layers) and H < H, (trajectory where the interaction between the 
motion with respect to the center of mass and the center of mass trajectory is taken into 
account). Hy, is the limiting altitude that is calculated to be 70-80 km.. In the region, 
H <= H,, the refined asymptotic method breaks down. An approximate solution is 
derived by formulating averaged equations for the asymptotic solutions that describe the 
change of the slow varying components of : 

1. The maximum and minimum values of the angle of attack during each 
oscillation 


2. The instantaneous oscillation period 


3. The averaged components of the center of mass trajectory 


The derivation of this method is mathematically rigorous. [Ref. 56:pp. 319-336] 
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The direct solution to this problem can be obtained by numerically 


integrating the complete equations of motion given by [Ref. 56:p. 45] 


d’a +My Gy Xypeeest) FER sNype-0) Se = 0 (125) 


dx; 
mr Bee (A et>5,--,0) = 0 








(126) 
(i = 1,2,...,6) 
where 
M, mG cr SI COS a) SCI) (127) 
SIN a 
X, = f(7,a)G + ((f(7,a) —f(7,@)] cosa -L(7,a) sina)r (128) 
pee f ATO) (129) 
X, = pcg eee sin@ (130) 
m 
Xx = 8, CcOsé _y (31) 
4 SP — 
V Rg, 
X, = -V sing (132) 
R 
X, = ae cos 6 (133) 
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ae: cosa -Q) Sin a 
I 


Zz 


G = \ sinta+r cosa = constant (134) 


However, in 1970, this numerical integration required an excessive amount of machine 
time or computer time, since the solution contains high frequency components along the 
trajectory. 
c. Follow-on Investigations of Uncontrolled Reentry Body Motion 

References [59] and [60] are extensions of Duzmak’s investigation. Each 
paper examines a certain aspect of the uncontrolled three-dimensional motion of a reentry 
body relative to its center of mass. 

Reference [59] investigates the three-dimensional uncontrolled motion of 
a spherical body relative to its center of mass with an arbitrary angle of attack. The 
mean differential equations are derived for rolling motion. These equations are solved 
in implicit form for any angle of attack by using the Van der Pol method of integration. 

Reference [60] investigates the three-dimensional uncontrolled motion 
relative to the center of mass of a reentry body with a small geometric and dynamic 
asymmetry. An approximate analytical solution for the equations for unperturbed three- 
dimensional and the averaged equations for perturbed three-dimensional motion are 
derived. By imposing no limits on the initial angular velocity and the three-dimensional 
angle of attack or nutation angle, the averaged equations of motion computational time 
is reduced by a factor of approximately three as compared to the numerical integration 


of the complete equations of motion. 
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B. SIX-DEGREE-OF-FREEDOM SIMULATIONS 

The foundation for six-degree-of-freedom motion has been discussed in the previous 
section. Several other investigators have examined the effect of motion relative to the 
center of mass on the trajectory of the center of mass [Refs. 61-64]. The equations of 
orbital decay, (121) through (124), when coupled with the equations of attitude motion, 
(112) through (114), will completely describe the motion of a vehicle during the final 
stages of life. Equations (125) and (126) are the coupled equations of six-degree-of- 
freedom motion written as a function of angle of attack. 

Coupling these equations means: the solution of one set of equations determines the 
magnitude and form of the forces or moments in the other set [Ref. 62:p. 13]. For 
example, if the area changes due to either a loss in mass or a change in attitude, then the 
ballistic coefficient (W/C,A) and the coefficient of drag will be affected. 

References [61], [62] and [63] investigate Skylab’s attitude and trajectory motion. 
Specifically, references [62] and [63] derive a variation of the coupled equations 
presented in the previous section. The authors performed a six-degree-of-freedom 
trajectory simulation using the coupled equations of motion. Figure 28 shows the flow 
chart used in their simulation [Ref. 62:p. 14]. Numerical integration of the six 
differential equations yielded altitude and orbital elements as a function of time. 
Computer run results for these simulations were extremely long and therefore the 
decision was made to artificially increase the magnitude of the aerodynamic damping. 
Subsequently, the results do not simulate the actual dynamical behavior, but they show 


the "possible" types of dynamical behavior for Skylab. [Ref. 62:p. 15] 
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Figure 28: Six-Degree-Of-Freedom Simulation Flowchart 
[Ref. 62] 
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Reference [64] is a reentry analysis of a proposed radioisotopic thermoelectric 
generator (RTG) that was connected to a generic user Satellite. A 3-degree-of-freedom 
reentry trajectory simulation was conducted that utilized aerodynamic, material property 
and heating characteristics. 

Newtonian and free molecular drag coefficients were calculated using the Mark IV 
Supersonic-Hypersonic Arbitrary-Body program for the generic satellite and each of its 
subelements. Figure 29 shows the three-dimensional shapes of the generic spacecraft and 
one subelement (reactor with radiator) used by the program [Ref. 64:p. 6]. Aerodynamic 
blockage was neglected and a zero angle of attack was assumed in the calculations. The 
heat transfer methodology used in the simulation to predict the satellite’s breakup during 
reentry, implemented a calibrated heat transfer model to closer simulate actual conditions 
{Ref. 64:pp. 7-12]. This will be discussed further in the section on breakup later in the 
chapter. The trajectory simulation was designed to change the mass properties and 
aerodynamic coefficients as the shape changed due to an "assumed predetermined 
breakup" scenario. This assumption was based on the basic geometric components or 
subelements of the satellite and their probable separation sequence during breakup. For 
example, the boom separated from the main part of the spacecraft and the reactor failed 
first, this was followed by the heat radiator cone, and then the other components in 
sequence. [Ref. 64:p.14] 

Trajectory simulations were run to assess the breakup altitudes, due to reentry 
heating, and minimum footprint lengths, due to fuel pin release at various altitudes, 


which were independent of the heating effects. 
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Figure 29: Simulated Generic Spacecraft And Reactor Sub-Element 
[Ref. 64] 
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C. RAREFIED GAS DYNAMICS 

Numerous authors have investigated reentry vehicle attitude, heating rates and 
aerodynamic coefficients in relation to the atmospheric flow regimes [Refs. 22, 56, 65- 
75]. As stated in Chapter II, the changes in flow regime and corresponding changes in 
critical parameters of the reentry trajectory or heating equations, are poorly understood 
in the transition from free-molecular to hypersonic continuum flow. For engineering 
applications the quantities of lift, drag and heat transfer are usually estimated by a 
"judicious faring" between regimes [Ref. 22:p. 203]. 

One solution to this problem is presented by Koppenwallner and Johannsmeier, 
reference [76]. This solution is a technique of "bridging" between the free molecular and 
continuum flow, based on Newtonian and free molecular theory. This technique allows 
the prediction of lift and drag during the hypersonic entry phase. [Ref. 76:p. 461] 

Three hypersonic flow regimes, with approximate boundaries, are described as 


follows: [Ref. 76:p. 461] 


1. Free molecular flow Kn > § 
2.  Rarefied transitional flow S > Kn > 0.001 
3. Hypersonic continuum flow Kn < 0.001 


These boundaries are actually dependent upon reentry vehicle shape and on the 
aerodynamic property considered. Figure 30 shows these flow regimes in an altitude- 


velocity graph [Ref. 76:p. 461]. 
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The analysis of this technique is limited to flow conditions where the hypersonic 
Mach independence principle is applicable; therefore, blunt shaped reentry vehicles will 
have a lower limit of Mach 5. Figure 31 shows the typical aerodynamic data variation, 
for a blunt body, through the three flow regimes of interest [Ref. 76:p. 461]. 

The typical drag coefficient behaves such that in free molecular flow, the drag 
coefficient is independent of Knudsen number and the value of Cp = 2. The drag 
coefficient decreases throughout the transitional flow regime and again reaches a constant 
value in the continuum flow regime. [Ref. 76:p. 461] 

The typical lift coefficient behaves such that in the free molecular flow, the lift is 
very small or negligible. The lift coefficient increases throughout the transitional flow 
regime until it reaches a hypersonic continuum flow value. [Ref. 76:p. 461] 

The heat transfer Stanton number (St) 1s very close to a value of one in the free 
molecular flow regime. In the transitional flow, the heat transfer efficiency decreases 
and the Stanton number decreases. In the continuum flow, the Stanton number decreases 
continuously as a function of Reynolds number in the manner St~ 1/VRe. [Ref. 76:p. 
461) 

Table III describes the bridging dependencies as modeled above. The derived 
transitional functions D, and L, are functions of angle of attack, vehicle shape and 
Knudsen number. These transitional functions must bridge the free molecular and 
continuum flow regimes and must degenerate, in the two limits, to the free and 


continuum flow values as shown in Table IV. [Ref. 76] 
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Figure 31: Aerodynamic Coefficients In The Flow Regimes 
[Ref. 76] 
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Table Ill: BRIDGING DEPENDENCIES 


[Ref. 76] 
a Sem 


a _ =— 


FLOW REGIME | DRAG | LIFT 
Cp = Ds (a, shape! 
Cp = D; (a,shape, Kn) _| C, = Ls (a,shape,Kn) 


Hypersonic Cy, = D, (a, shape) C, = L, (a, shape) 
continuum 





Table IV: TRANSITIONAL BRIDGING FUNCTIONS 
[Ref. 76] 


) DEGENERATE TO TRANSITIONAL DEGENERATE TO 
FREE MOLECULAR BRIDGING FUNCTION | HYPERSONIC CONT 


D, = (a, shape) D, = (a, shape, Kn) D, = (a, shape) 
L, = 0 


e 
s Ud 


L, = (a, shape, Kn) L, = (a, shape) 





Local flow independence exists in free molecular as well as hypersonic Newtonian 
flow. This means that unless "shadowing" exists, surface and shape elements will not 
influence each other. The shape elements of this technique are cones, spherical caps and 
cylindrical shells. These shape elements allow the composition of capsules, blunted 
cones and cone-cylinders. Figure 32 shows the shape elements and several composed 


bodies [Ref. 76:p. 462]. 
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Figure 32: Shape Elements And Composed Bodies 
[Ref. 76] 


116 


1. Newtonian Aerodynamics for Hypersonic Continuum Flow 
For the case of blunt bodies, it can be assumed that the main contribution to 
drag is pressure; therefore, the influence of viscous effects on the aerodynamic 
coefficients are neglected. [Ref. 76:p. 462] Using the Newtonian pressure law, the 
coefficients of drag and lift may be described in differential equations. The basic 


Newtonian relationships are as shown in Table V. [Ref. 76:p. 462] 


Table V: NEWTONIAN LOCAL PRESSURE LAW 
(Ref. 76] 


[Wetted surface [lal < 90° pda = ky c0s%0 
[Newtonian wake | lel > 90: p | On = Coy 








where 
9) = pressure 
Go» = dynamic pressure in free stream 
k, |= Newton factor 
6 = flow inclination against surface normal 
K = ratio of specific heats 
Cp. = coefficient of drag as a function of angle of attack 
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and 
ky =2 (simple Newton) 
k, = (x + 3)/(x + 1) (modified Newton) 


In the Pike formulation, the Newtonian differential equation is given by [Ref. 71:p. 462] 


A, (a) 


Cy (a) + Cp (a)cot(a) + IZC Oe 0K. (135) 





R 
3C,(a) = Cy (136) 


where 
A, = Flow projected area of body 
A, = Reference area (xD*/4) 
In order to determine the drag coefficient and solve the differential equations, 
the following approach is taken: Cy) (a=0°, shape). This implies that the zero angle 
of attack drag is proportional to the Newton factor and to a shape dependent factor, C, 


Ref. 76:p. 462] 
C,(a=0°) = k, C; (137) 


Upon integration of the Newtonian pressure distribution, the following values 


for C, are determined for various bodies in Table VI. [Ref. 76:p. 462] 


118 


Table VI: SHAPE FACTOR VALUES 
[Ref. 76] 


= —————E= ——— 


a Cr 

Sphere 

Cylinder p= 0.67 Ieross flow) | 
mea 






Spherical cap Cys= 1-1/8 + (d/ty)? 


Figure 33 shows the shape factor as a function of the ratio (d/r,), diameter 


over nose radius [Ref. 76:p. 466]. 

When the shape factor or drag at zero angle of attack are known, the 
aerodynamic coefficients of the class are completely fixed. The following universal 
solutions for drag and lift are obtained from the Newtonian differential equations [Ref. 


76:p. 462) 


119 


1.0 
: an i 
? Cs s1e1.(d) | 


07 ‘ 
7 toh cal am 
OS 


Figure 33: Shape Factor As A Function Of (d/r,) 
[Ref. 76] 
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k * 

C= (2c, - (5C, - 3) sin%(a)} cosa (138) 
K e ° 

7+ ( 2 (1 - 2C,) - (5C, - 3) sin*(a)) sina (139) 


These solutions are valid fora < a,,, , since at a,,, the wetted area decreases due to 
Newtonian shadowing. Figure 34 shows the universal Newtonian lift and drag functions 
for various angles of attack [Ref. 76:p. 466]. 

The shape factor at zero angle of attack serves as the critical parameter in this 
method. The solution is valid only under the condition that the wetted surface area 
remains constant for any angle of attack. [Ref. 76:p. 462] 

As the geometric body bluntness increases, C, increases and Newtonian 
shadowing is shifted to higher angles of attack. Table VII shows that depending upon 
the body shape, the lift slope at zero angle of attack may be either positive or negative. 


[Ref. 76:p. 462] 
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Figure 34: Newtonian Lift And Drag Functions 
[Ref. 76] 
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Table VII: BODY SHAPE AND LIFT SLOPE 
[Ref. 76] 


Body shape Lift slope @ a=0 





This demonstrates that blunt reentry bodies will experience a negative lift for a positive 
defined angle of attack. [Ref. 76:p. 462] 
2. Free Molecular Flow Model 

The technique of reference [76] also develops analytical formulas for free 


molecular flow. The general considerations include: [Ref. 76:p. 463] 


1. Accommodation coefficient, o 
2. Finite molecular speed ratio, S 


3. Wall temperature, Ty 


These factors alone are insufficient thus the following simplifying assumptions are made: 


[Ref. 76:p. 463] 


1. Body shape Hypersonic blunt, Ma (d/l) > 2 
2. Aerodynamic cold wall Tae =< J 
3. Diffuse molecular reflection og = 1 
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From these constraints the local surface pressure and shear stress may be 
defined as shown in Table VIII. The aerodynamic force coefficients are also derived as 
Shown in Table IX. 


Table Vill: LOCAL SURFACE PRESSURE AND SHEAR STRESS 


[Ref. 76] 
a 


lal < 90° pid. = 2 cos%6 
lal < 90° rig..= sind-cosd 
















Table IX: AERODYNAMIC FORCE COEFFICIENTS 


(Ref. 76] 


Drag coefficient Cy= 2°A,(a) / Ap(O) = cos a 
Lift coefficient 









These assumptions imply that no lift is produced and drag is proportional to the flow 
projected area, Ap(a), in the free molecular flow regime. [Ref. 76:p. 463] 

In the more general free molecular flow case, where T,,/T,, , o, and o, are 
considered, an axisymetric blunt body formulation for the force coefficients may be given 


by 
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x 17, | ] 
C, = 20,cosa + a, . Tr [cues | VE Eaed cost - (140) 
(2 -0,-0) =[3(C,.+1)cose + (5Cy_~3)cos3ar 
yx | 7 
nx | 4w 4 
CG, oS Tt c, 3] sin2a + (141) 
(20: -0)- | ~(C,,+1)sina -(5C,,-3)sin 30 
where 

T,. = temperature, free stream 

Cs, = shape coefficient, in front of normal shock wave 

C,, = shape coefficient, behind normal shock wave 


The physical significance of several of the terms in equations (140) and (141) 

are: [Ref. 76:p. 463] 
1. The first term in Cy gives, for diffuse reflection, o=1, the contribution of the 
incident flux to the aerodynamic coefficients. 


2. The T\/T, term states the influence of the reentry vehicle surface wall 
temperature on drag and lift. 


3. The (2-o,-0,) term vanishes for diffuse molecular reflection. 


In simple free molecular flow the equations will degenerate appropriately as 
6, =1, o,=1 and T,/T,=0. In the Newtonian formulation the equations will degenerate 
as o,=1, o,=0 and T,,/T,.=0. In this last case C,, vanishes and C,, is the Newtonian 


shape coefficient, C,. [Ref. 76:p. 463] 
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3. Bridging Free Molecular and Continuum Flow 

Several methods exist which attempt to bridge the gap between the free 
molecular and continuum flow regimes. [Refs. 77-79]. In the former USSR and DLR, 
local bridging is accomplished with a finite surface element method, [Refs. 77-78]. In 
the U.S., reference [79], bridging has been accomplished through an integral coefficients 
method. [Ref. 76:p. 463] 

The method of bridging by shape element description is developed in 
reference [76] and is presented in the following section. The methods of references [78] 
and [79] may be used to derive analytical formulas for trajectory calculation. However, 


the basic bridging relations must be derived through experimentation. [Ref. 76:p. 463] 


a. Shape Element Bridging Method. 

Experimental data, presented in Figures 35 and 36 show the drag 
coefficient changes of a sphere and disk respectively. The data covers the the entire 
transitional flow regime [Ref. 76]. 

For a sphere, the Reynolds number, behind the normal shock (Re,), can 


be related to the Knudsen number, in the free stream (Kn,), as follows [Ref. 76:p. 463] 


Res = 1.26) ee (142) 
k-1 | Ky 


From this experimental data, an approximate formula for Cp, as a 





function of Reynolds number and shape, can be derived by using the "reduced" 


aerodynamic coefficients. [Ref. 76:p. 463] 
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Figure 35: Sphere Drag Coefficient In Rarefied Flow 
[Ref. 76] 
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Figure 36: Disk Drag Coefficient In Rarefied Flow 


[Ref. 76] 
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C, = ColRé) “Foe (143) 
Corm~ pc 
UD = L/D(a=const., Re.) (144) 
L/D(a@ =const., Continuum) 
where 
Coc = drag coefficient, continuum 


Corm = drag coefficient, free molecular 
Cp = reduced drag coefficient 


L/D 


I 


reduced lift-drag ratio 
Analysis of the reduced coefficients at zero angle of attack shows the 


following: [Ref. 76:p. 463] 


1. Slope of the reduced drag coefficient is shape independent. 
2. Continuum and free molecular flow boundaries are shape dependent. 


3. Continuum and free molecular drag coefficients are shape dependent. 


Therefore the bridging function is derived from the reduced coefficients as follows [Ref. 


76:p. 464] 


1 
. 3205 | Coru~ 


C, = (145) 





And from this, the drag coefficient may be extracted in one of two forms 


Ie, 


Cy = =a ne (Com ~ Coo) In = * es (146) 
c.g ON) a a (147) 
C PaGN ep a6) anee 


The shape dependent boundaries are defined as: 


1. Ree > Re, > Rey 


2.” Knee Kn Skee 


and the experimentally derived values are summarized in Table X below. [Ref. 76:p. 


464) 


Table X: SHAPE DEPENDENT BOUNDARIES 


(Ref. 76] 


EE Continuum Flow Free Molecular Flow 
Se Re@5¢ Kn Rory Kren | 


Disk Disk § ~— | 2090 0.117 | 0.115 2185 


Sphere 2 | 89.3 0.037 | 0.492 6.8 





In summary, the technique of reference [76] allows the following 


conclusions to be drawn: [Ref. 76:p. 464] 


1. Newtonian theory based on Pike’s method for the hypersonic continuum flow 
is useful, as it shows the shape dependent and shape independent aerodynamic 
coefficients. 


130 


2. The bridging function is based on experimentally derived results for the 
transition flow. 


3. Lift and drag bridging do not follow the same laws. 


b. Local Bridging Method 
Reference [78] develops a technique of "local" bridging. The differences 
between the shape element method, referred to as a “global” method, and local bridging 
are: [Ref. 78:p. 469] 
1. Global bridging performs across the spectrum of aerodynamic coefficients for 
the complete body. 
2. Global bridging is usually tailored to a specific class of shapes. 


3. Global bridging techniques require new experimentally derived "fitting 
constants” for each new shape of interest. 


4. Local bridging is a method of bridging the local pressure and shear stress 
coefficients. The local distribution is then integrated over the body surface 
which yields the global force and moment coefficients. 

The common principle of all bridging techniques is the manner in which they model the 
transition function. The known free molecular and continuum flow limiting values for 


a specific aerodynamic coefficient (lift, drag or heat transfer) are weighted and applied 


to the bridging function. The general case is as follows 


C(X) = Cry AX) + Co (1 -f)) (148) 
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where 


@ 
| 


aerodynamic coefficient considered 


~ 
I 


rarefaction-dependent flow parameter 
The potential improvements of local bridging over global bridging are: 


[Ref. 78:p. 470] 


1. Ease of adaptation to different shapes without changing internal constants. 
2. More reliable calculation of moments. These are very sensitive to rarefaction 
because of the varying contributions of pressure and shear in the different flow 


regimes. 


3. Simplified method for determining reference quantities for local coefficients, 
such as reference length based on a local coordinate. 


4. Ability to account for non-uniform flow field conditions. 


The conclusion of reference [78], after examining several different 
bridging methods and various shapes, is that the current state of experience in applying 
the local bridging methods allows no definitive answer as to whether or not they can 


serve as an effective analysis tool in transition flow analysis. 


4. Gas-Surface/Gas-Gas Interactions 
Gas-surface interactions are significant in the understanding of reentry 
dynamics and aerothermodynamics. The primary cause of concern is that at orbital 
altitudes, the highly rarefied flowfield is dominated by gas-surface interactions that occur 
at average velocities corresponding to that of the orbiting vehicle. [Ref. 80:p. 1] Under 


these reentry conditions, gas-gas interactions become important and gas molecules that 
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reach the vehicle surface tend to have lost some of their initial translational energy. This 
loss of translational molecular energy 1s due to conversion into other forms such as heat, 
internal energy, etc. Figure 37 shows molecular velocity distributions in rarefied and 
transitional flow about a reentering sphere [Ref. 80]. The collision mechanisms by 


which the translational energy is converted include: 


1. Chemical reactions 
2. Ionization 


3. Dissociation 


Since the nature of the gas-surface interaction is known to be dependent upon the velocity 
and energy of the incident molecules, it becomes necessary to know the state of the gas 
molecules reaching the surface. [Ref. 80:p. 1] 

Wilmoth, et al, [80] uses the technique of Direct Simulation Monte Carlo 
(DSMC) as developed by Bird, [Ref. 81] where the molecular velocity and energy 
distributions of the gas molecules are a direct result of the simulation process. Bird uses 
a simple engineering model of the gas-surface interactions, which accounts for diffuse 
and specular reflection along with other phenomena. This model can accommodate 
processes such as catalytic reactions and molecular recombination. A limitation of this 
method is that parametric studies must often be performed in order to place bounds on 
the predicted quantities of interest. [Ref. 80:p. 1] This also implies that the analyst must 


judiciously apply the model based on experience and a limited base of experimental data. 
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Figure 37: Molecular Velocity Distribution 
[Ref. 80] 
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The lack of experimental data is the reason why more detailed gas-surface 
interaction models have not been developed. This results mainly for two reasons: [Ref. 
80:p. 1] 

1. It is very difficult to simulate gas-surface interactions at orbital or entry 
velocities in a laboratory. 

2. It is difficult to characterize the surfaces used in laboratory experiments with 
sufficient generality that the results may have application in an engineering 
context. 

Because of the tendency of gas molecules to decelerate after gas-gas collisions, which 
occur before reaching the vehicle surface, it is important to quantify the actual velocities 
and energies encountered in such cases. 

Reference [80] studies the "typical" gas-surface interactions for transitional 
flow at entry velocities. The reentry vehicle was a 1.6 m diameter sphere at a free 
stream velocity of 7.5 km/sec over an altitude range of 90 to 130 km. [Ref. 80:p. 2] 
The study was conducted using the DSMC method of Bird. 

The flow conditions of the DSMC simulation are given in Table XI [Ref. 
80:p. 6]. The atmosphere is modeled by Jacchia, 1977, with an exospheric temperature 
of 1200K. The surface temperature of the sphere was assumed constant at 350K. The 
gas-surface interaction was assumed to be diffuse with full thermal accommodation, and 
the surface was non-catalytic. Five atmospheric molecular species were modeled O,, N,, 


O, N and NO, with 23 reaction possibilities. [Ref. 80:p. 3] 
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Table XI: FREESTREAM CONDITIONS 
(Ref. 80] 


=e EP EEE 
kg/m° km/s Fraction He A ins g/mol m 
[0.700 
BoE Se 
Fro | senor | rs | zer | otas | 0.70 | 0.100 | 27.22 | 0.500 
p10 | zanior_| 7s | 202 | ooss_| 0.709 | o.199 | 20.14 | 2.001 
130 | ezaor | zs | soo [oom | oser [ozs | 25.00 | 7.724 




























Wilmoth’s results show the following: [Ref. 80:pp. 3-5] 


1. At 130 km, a density rise of nearly 22 times the freestream value occurs over 
a distance of ~3 m. However, based on analysis of the data it is determined 
that very few collisions are occurring. 

2. At90km, the density increases to well over 100 times the freestream value over 
a distance of ~0.1 m. In this altitude range collisions become very significant 
and there is considerable chemical activity. 

3. There are significant variations in the velocity and energy of molecules reaching 
the surface over an altitude range of 130 to 90 km. Figure 38 shows the 


average translational energy per particle striking the reentry vehicle surface 
[Ref. 80]. 


D. SURFACE ROUGHNESS EFFECTS 

Another important consideration in the reentry phase is the "boundary-layer 
transition.” It is this phenomenon which is responsible for heat transfer to/from the 
reentry body due to atmospheric contact with the body. In the continuum flow regime 


viscous effects are essentially restricted to a small layer called the boundary layer. It is 
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within this boundary layer that the details of "fluid" motion determine the levels of skin 
friction and heat transfer from the flow. [Ref. 22:p. 211] 

A parameter commonly used for characterizing the boundary layer is the Reynolds 
number. [Refs. 82-83] The Reynolds number is defined as the ratio of inertial forces to 


the viscous forces, as given by [Ref. 22:p. 212] 


inertia forces _ A(momentum)/(unit time) _ Re (149) 
viscous forces (shear stress)x(unit area) 


In order for a fluid (the atmosphere 1s modeled as a fluid in the continuum flow) 
to support a shear stress there must be relative motion between adjacent layers. This 
implies that there is a velocity differential within the flow layers. The following terms 


are defined: fRei..22:n oy 


T = shear stress = p(dV/dy) 
fe = dynamic viscosity 
y = coordinate normal to direction of motion 


The definition of Reynolds number may now be given as [Ref. 22:p. 213] 


_ dmvyidt _ [pL3v/(Liv)] _ jw. ML (150) 
u(dV/dy)A u(V/L)L? fr y 
where 
yp = kinematic viscosity = (p/p) 


Now the Reynolds number may be used to characterize the boundary layer flow 
conditions. When the viscous forces are large enough to damp out the oscillations caused 


by the dynamic forces, the flow is laminar and the Reynolds number is small. 
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Conversely, the flow is turbulent when the dynamic forces overcome the viscous forces. 
Also, the Reynolds number is large for turbulent flow. [Ref. 22:p. 212] 

The velocity profiles within the laminar and turbulent boundary layers show some 
significant differences. The magnitude of velocity in the turbulent boundary layer is 
notably greater, especially near the reentry body’s surface. This implies that the 
turbulent boundary layer yields much more energy near the surface than does a 
corresponding laminar flow. The velocity profile determines the skin friction on the 
surface and it can be expected that greater heat transfer will occur under conditions of 
a turbulent boundary layer [Ref. 22:p. 213]. Figure 39 shows the boundary layer types. 
Figure 40 shows boundary layer velocity-distance profiles. Figure 41 shows the altitude, 
air speed and dynamic pressure [Ref. 22:pp. 213-214]. 

A Reynolds stress turbulent boundary layer model which specifically accounts for 
surface roughness effects is described by reference [84]. In this study, surface roughness 
is represented by distributed "sources" and "sinks" in the various governing equations. 
The most significant term is a sink term in the mean momentum equation, which 
represents "form drag" on the roughness elements. [Ref. 84:p. 2] 

A fundamental assumption of this model is that the flow around the individual 
roughness elements (only distributed roughness 1s considered) is attached to the elements. 
[Ref. 84:p. 3] The roughness elements provide a distributed sink, due to drag, for the 
momentum equation and distributed sources for mean turbulent kinetic energy and 
dissipation. This model also assumes that the roughness elements occupy no volume; 


therefore, this assumption becomes more severe as the roughness density increases. In 
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order to compensate for this, the model has been extended to account for the blockage 
effect of the roughness elements. [Ref. 84:p. 4] 


The mean momentum equation is given as [Ref. 84:pp. 4-5] 


fo) pUL +pv 5 = - fy Pr 2 | 2 
Ox Ox oy Poy (151) 
——, 1 > D aD? )~ 
sai =p C= | Pe 
50 Ae oe | 


The variables U, V, u/ and v/ are the reduced variables, under the boundary layer 
approximation, as given in reference [85]. These variables and this equation will not be 
derived here for the sake of brevity. 


The major advantages of this model are: [Ref. 84:p. 5] 


1. Solutions are obtained for both velocity and thermal variables. 

2. Heat transfer is obtained directly, without invoking Reynolds analogy. 
3. Finite difference solutions are obtained using the boundary conditions that 

fluctuating quantities are zero at the solid wall and in the free stream. 
Reference [84] concludes with a comparison of a smooth wall turbulent boundary 
layer model and the developed rough wall model. The rough wall model is determined 
to show that roughness spacing is more critical than roughness height, under the 
conditions tested; however, the limited skin friction data obtained in the study cannot be 


interpreted unambiguously. [Ref. 84:p. 38] 
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E. REENTRY HEATING EQUATIONS 
In Chapter II, the fundamentals of reentry heating were discussed as presented by 
the authors of references [20], [40] and [43]. Aerodynamic heating, as it applies to the 
reentry of space vehicles, takes its roots in the work of Allen and Eggers, reference [65]. 
Reentry heating becomes an important consideration of the overall reentry process for 
the following reasons: [Ref. 42:p. 139] 
1. Structural performance of the reentry vehicle is dependent upon the dynamic 
pressures encountered during the reentry, which is a function of the reentry 
trajectory. 


2. Structural strength of the vehicle is a function of the stresses induced by 
temperature gradients within the component materials. 


3. Temperature gradients aré proportional to the time rate of heat input and 
maximum time rate of heat input. 
Therefore, the three critical parameters of the reentry trajectory are the total heat input 
along the trajectory, the maximum rate of aerodynamic heating and the maximum 
dynamic pressure. [Ref. 42:p. 139] 

The mechanism of heat flow into the reentry vehicle during atmospheric entry was 
first described in reference [65]. Since then, numerous combinations of reentry speed 
regimes and aerodynamic shapes have led to the publication of numerous technical 
reports. However, the basic aspects of aerodynamic heating during reentry are still the 
same. The numerical factors for different heat transfer formulas and their ranges of 
validity in terms of the regime of speed are the only variation among the numerous 


authors. [Ref. 42:p. 139] 


144 


The basic equations of reentry heating, equations (18), (19) and (20), are developed 
with their simplifying assumptions as presented in Chapter II, pp. 46-48. These 


assumptions yield the following limitations: [Ref. 42:pp. 141-142] 


1. Neglecting radiative heat transfer from the vehicle or to the vehicle from the 
high-temperature air between the shock wave and the vehicle surface, is based 
on the fact that the maximum allowable surface temperature is about the same 
for a variety of reentry vehicle shapes. Thus, outward radiation from the 
surface will be about the same. Neglecting the radiative heat transfer from the 
disturbed air is a qualitative simplification and therefore negates the application 
of the equations to very blunt and heavy shapes at reentry speeds of 3 km/sec 
or greater. 


2. Neglecting the real-gas effects in the flow, most importantly dissociation, on 
convective heat transfer is a good approximation for reentry speeds of up to 3 
km/sec. Nevertheless, this assumption is conservative and results in higher 
calculated heating rates than actual rates. 

3. Neglecting the shock-wave boundary-layer interactions implies that the laminar 
skin-friction coefficient, on a flat plate at zero incidence, is being held constant. 
This assumption is not valid at reentry speeds over 6 km/sec. 

4. Assuming a Reynolds analogy and holding the Prandtl number constant also 
restricts the validity of the equations to reentry speeds of less than 3 km/sec. 

For the case of low earth orbits and naturally decaying satellites, these assumptions 


create severe restrictions. The circular orbital velocity may be calculated from [Ref. 


Ds OO] 


v= |# (152) 
r 
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where 


ie radius from center of Earth 


V. circular velocity 


Using this equation, it 1s possible to approximate the circular orbital radius of a 
satellite traveling at an altitude of near 120 km or radius of 6499 km. This is the altitude 
of concern for the focus of this thesis and it can be shown that the orbital velocity of 
interest is approximately 7.8 km/sec. The orbital radius of a satellite travelling at speeds 
of 3 km/sec and 6 km/sec are 44,289 km and 11,072 km respectively. This equates to 
circular orbital altitudes of 37,910 km and 4,694 km . 

The assumptions of reference [65] which are maintained in references [14], [40] 
and [43] should be removed for an accurate quantitative analysis of the aerodynamic 


heating during the reentry of a specific vehicle. [Ref. 42:p. 142] 


F. STRUCTURAL BREAKUP OF A REENTRY BODY 

The previous sections of this chapter have served to develop the "state-of-the-art," 
as determined through the literature survey, of reentry formulations and solutions. The 
culmination of the uncontrolled reentry process is usually the structural breakup of the 
reentry body. The breakup of a reentry body is stated to be a function of surface 
temperature, in that, structural failure (breakup) is assumed or expected to occur when 
the outer structure reaches its melting temperature. [Refs. 63, 86] It has been shown 
previously that the maximum temperature a reentry body will experience is determined 


by the maximum heating rate, which is a function of ballistic coefficient. [Ref. 86:p. C- 
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11] Although the peak heating rate increases as the ballistic coefficient increases, 
average heating rate (a measure of overall survivability) and maximum local stagnation 
region heating (a measure of local “hot spots"), must be considered in determining where 
breakup occurs. [Ref. 40:pp. 181-182] 

In the 1970’s, the Air Force conducted reentry experiments that used optical and 
radar techniques to observe actual breakup events. The objective of these tests was to 
determine the survivability of reentry body debris. Specific findings from these tests 
indicate: [Ref. 86:pp. C-4--C-10] 

1. Classical convective heat transfer analysis underestimates the reentry body 
survivability. Specifically: actual breakup is at an altitude of at least 10 nm 
lower than predicted; actual surface temperature is at an altitude 10 nm lower 
than predicted; and the effective heating rate input is a factor of four lower than 
predicted. 

2. Consistent catastrophic failure of magnesium/aluminum structures is at an 
altitude of 42 nm. Figure 42 [Ref. 86] shows convergence of three ballistic 
coefficient lines close to the magnesium/aluminum melt zone at this approximate 


altitude. 


3. Phenomenon of breakup process is independent of: body attitude and rates; body 
diameter; body shape; and entry flight path angle (0.3° > 1.5°). 


4. Ballistic coefficient and material of a body determines survivability. 


5. A body with low ballistic coefficient will survive reentry, as will a body with 
a higher melting temperature survive reentry with a higher ballistic coefficient. 


6. Reentry body structural integrity is maintained until melting temperature is 
achieved. 


7. Surface structure temperature is determined by radiation equilibrium that is 
based on a shallow path angle and a low thermal capacity of the outer structure. 
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Figure 42: Structural Breakup - Heating Rate vs Altitude 
[Ref. 86] 
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Reference [85] concludes that classical convective heat transfer analysis in a free- 
molecular flow regime is not indicative of the transition/continuum flow heating which 
is responsible for structural breakup. This finding is in agreement with reference [42] 
which, as shown previously, discounts the classical reentry heating equations based on 
the limitations imposed by the simplifying assumptions. [Ref. 42:pp. 141-142] Reference 
[64] also cites the 1970’s Air Force studies and indicates the deficiencies of the classical 
convective heat transfer analysis. 

The consistent underestimation of survivability was further investigated in reference 
[64]. By using a calibrated heating model, the observed breakup altitude of 
approximately 42 nm (79.5 km) was successfully predicted for satellites with aluminum 
structures. Figure 43 [Ref. 64:p. 48] shows the predicted breakup altitudes for two 
satellite reentry simulations. [Ref. 64:pp. 7-13] 

The heating model was interfaced with a trajectory simulation model in order to 
estimate the altitudes where possible breakup events could occur. Table XII shows the 
breakup analysis results of the critical elements with their associated material 
composition, heating rates and breakup altitudes [Ref. 64:p.59]. 

Since the breakup analysis was assumed with strictly aerodynamic heating for 
structures with low melting points, the authors recommend further investigation for the 
combined effect of both aerodynamic heating and aerodynamic loading on breakup. 
Specifically, for the case where higher melting point materials will survive to altitudes 


where deceleration loads are significant. [Ref. 64:p. 21] 
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Figure 43: Breakup Altitudes For Corrected Heating Rates 
[Ref. 64] 
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Table XIl: RESULTS OF BREAKUP ANALYSIS 
[Ref. 64] 


onfiguratioa Critical Candidate Greq'd Q Breakup 


Construction Altitude 
Material(s) BTU/ Ec@=s) (3TU/ tees) 


Radiator 


Reactor/Radiation| Radiation 
Shield Shield 


Fuel Pin #1 115 = 1410 
vel Pin #2 247 = 1185 





a Breakup anticipated for tumbling mode, some 
pins may survive. 


b —s— Breakup anticipated in tumbling mode. 
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Another series of investigations of the breakup process are references [10] and [63]. 
In these studies, the reentry/breakup of Skylab was reconstructed by piecing together all 
available data - after the fact. Both of these ae show conclusively, through telemetry 
analysis, that the survivability of Skylab was underestimated. Reference [10] cites initial 
breakup predictions of 120 km and shows that breakup did not occur until at least 100 
km. References [63] and [10] disagree on exactly where the OWS SAS (Orbital 
Workshop Solar Array System) separated from the main body, based on telemetry 
received at Ascension Island. Reference [10] states that the array was intact over 
Ascension Island; however, it was not generating its predicted output. Therefore, it was 
concluded that the array was either bent back or physically damaged. 

Reference [63] concludes that the OWS SAS was completely separated from the 
main body prior to telemetry acquisition at Ascension Island [Ref. 63:p. 344]. Finally, 


reference [63] postulates a probable breakup scenario as follows: [Ref. 63:p. 344] 


1. OWS SAS array (aerodynamically) off at 62 nm / 117 km. 
2. ATM separates from remaining OWS at 54 nm / 102 km. 
3. ATM SAS arrays separate between 50 and 54 nm / 94.5 and 102 km. 


4. ATM and OWS breakup at 42 nm / 77.8 km. 


Reference [10] concludes that breakup did not start until some altitude below 100 km. 
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IV. DETERMINISTIC REENTRY/IMPACT PREDICTION METHODS 


A. CURRENT REENTRY/IMPACT PREDICTION METHOD 

The following section of this chapter describes the various techniques or methods 
used by different countries or organizations dealing with reentry and impact prediction 
of naturally decaying objects. It is therefore useful to the reader to understand the 
current U.S. method used in reentry/impact prediction, since it is the standard by which 
all other methods (as determined by this literature survey) are compared. 

As Stated in previous chapters, the current method for predicting reentry time and 
impact location is that of the Space Surveillance Center (SSC) located at Cheyenne 
Mountain Air Force Base, or commonly referred to in the literature as NORAD. [Ref. 
87] NORAD produces "element sets” which are mean values of the orbital elements that 
have been obtained by removing the periodic orbital variations in a particular manner. 
In order to use these element sets, and obtain reasonable predictions, these periodic 
variations must be "reconstructed" by the prediction model in precisely the same manner 
as they were removed by NORAD. Therefore, an input of NORAD element sets into 
another model (even though it may be more accurate, or even into a numenical integrator) 
will result in degraded predictions unless, as previously stated, the new model can 
"reconstruct" the periodic variations. [Ref. 87:p. 1] 

NORAD element sets are generated with a general perturbations (GP) model called 


SGP4. SGP4 was developed by Ken Cranford in 1970. This model was a result of 


WSS 


simplification of the more extensive analytical theory of Lane and Cranford (1969) which 
uses the solution of Brouwer (1959) for its gravitational model and a power density 
function for its atmosphere model. [Ref. 87:p. 3] It should be emphasized that this 
atmosphere "model" is a Static representation of density as opposed to the dynamical 
models discussed in Chapter II, such as those of Jacchia and others. 

The gravitational model includes J, and J, harmonics; however, J, and J; were 
dropped in order to avoid singularities occurring at critical inclination [Ref. 88:p. 2]. 
Rates of change of mean motion and eccentricity are derived from the density function. 
The product of ballistic coefficient and a reference density, denoted B*, is treated as a 
solved for parameter. Coupling between J, and drag is included in the argument of 
perigee, right ascension of the node and mean anomaly. The mean motion of SGP4 is 
a pure Brouwer, or two-body, mean motion. [Ref. 88:p. 3] 

Reentry and impact predictions, in the U.S., are made using a special perturbations 
(SP) propagator with conversions between GP and SP theories handled as outlined in 
references [87, 89-90]. The GP theory is fast, analytical and of low-accuracy, when 
compared to the SP theory. SP theory uses a Gauss-Jackson, eighth-order, numerical 
integrator, incorporates a 6-12th order geopotential model and applies a dynamical 
atmospheric density model (Jacchia-65). The "conversion" between GP and SP theories 
is the process of performing an SP differential correction of the initial state vector as 
derived from the GP theory (NORAD two-line element set). These are the initial 
standards for GP and SP theory compatibility which must be considered in the following 


discussion of different reentry/impact prediction methods. [Ref. 91:p. 8] 
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The first TIP (Tracking and Impact Prediction) run 1s performed approximately 15- 
20 days prior to the estimated reentry date, which is initially predicted by the GP model. 
This is also when tasking of the observation sites is initially increased in order to support 
high accuracy SP processing. [Ref. 92:p. 1] 

Orbit determination 1s accomplished through a first-order, linear, weighted, least 
Square, curve fitting process, commonly called differential correction. Sliding fits are 
used to process both new and old metric data until the satellite is "no longer in orbit." 
The force models used are the Jacchia dynamic atmosphere (1965) and the World 
Geodetic System Earth gravity model (1972). The Earth gravity model is truncated to 
the sixth order for satellite decay predictions. [Ref. 92:p. 2] 

Direct step-by-step numerical integration of the total acceleration acting on the 
decaying satellite is accomplished in the manner of Cowell’s method. Gauss-Jackson 
eighth order predictor-corrector formulas, in ordinate form, are used to integrate the 
equations of motion. Because of computer run-time constraints, the partial derivatives 
necessary for differential correction are computed analytically except for the secular 
variations due to atmospheric drag of the orbit semi-major axis and eccentricity. These 
parameters are integrated numerically using a low-order (trapezoidal rule) integrator. 


[Ref. 92:p. 2] 
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Currently, the differential correction solution state consists of the equinoctial 
elements and satellite ballistic coefficient model parameter, By. After each solution, the 


new state is used to predict a decay time (when altitude equals 10 km). [Ref. 92:p. 2] 


B. ALTERNATE REENTRY/IMPACT PREDICTION METHODS 


1. Reentry Prediction Methods At ESOC 

One method used at the European Space Operations Center (ESOC) for the 
prediction of reentry and impact of decaying satellites is an improved and computerized 
version of the King-Hele technique, reference [55]. [Ref. 93] Another method used at 
ESOC is based in a computer program called FOCUS [Ref. 94]. 

The principal characteristic of FOCUS is its ability to overcome the 
deficiencies of the semi-analytical orbit prediction techniques at low altitudes [Ref. 94:p. 
26]. Recall that in near-Earth orbits, 200-700 km, Earth oblateness (J,) is regarded as 
the only first-order perturbation. Higher zonal harmonics and air drag are regarded as 
second-order contributions. All other effects are considered as less than second-order 
perturbations (less than J,” in magnitude). [Ref. 94:p. 25] Eventually, the satellite passes 
through an altitude regime where the air drag force (caused by increasing atmospheric 
density) reaches a magnitude of the same order as the J, Earth oblateness effect. At this 
altitude, = 150 km, depending on the vehicle’s ballistic coefficient, the accuracy of 
analytically derived drag perturbation results strongly deteriorates. [Ref. 94:p. 26] 

FOCUS stops the state propagation after passing through a user-defined 


altitude shell (h = 170 km specifically for the case of Salyut-7/Cosmos-1686, which is 
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the focus of reference [94]) and forwards a calculated osculating Keplerian state vector 
at epoch, together with all relevant perturbation parameters. This new set of data is then 
numerically integrated and the reentry trajectory 1s propagated until shortly before impact 
with the Earth. [Ref. 94:p. 26] 


Several key features of the FOCUS program are: [Ref. 94:p. 26] 


1. Perturbation equations: 

(a) Cowell’s formulation of the perturbed Newtonian equations written in terms 
of six first-order, differential equations for each component of the Cartesian 
State vector 

(b) Reference frame is the mean equatorial system of date 

2. Perturbation models: 

(a) Geopotential model GEM 10B (J,-J, used) 

(b) Atmospheric density models: MSIS-86 for altitudes = 120 km, U.S. 
Standard Atmosphere (USSA-76) for altitudes < 90 km, and a bridging 
function for altitudes between 90 and 120 km. 

(c) Variable drag coefficient, Cp) = f(Ma, Re, Kn) 

(d) Co-rotating Earth atmosphere 

(f) Luni-solar third body attraction (point mass) 

(g) Solar radiation pressure 

3. Integrator: 


(a) Runge-Kutta/Shanks 7/8 single step method for generation of a starting arc 


(b) Adams-Bashforth/Adams-Moulton(AB/AM)forth-orderpredictor-corrector, 
multi-step method for propagation of the Cartesian state vector. 
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(c) Non-regularized time (t) used as an integration variable, with constant step 
sizes of At = 30 sec. 

A significant limiting feature of this program is that the reentry trajectory is 
terminated at 30 km altitude because the governing laws of perturbed Keplerian motion 
become invalid below this threshold altitude. This criterion is marked by a decrease in 
the orbital energy to a level, where the aerodynamic forces are essentially in balance with 
the zero-th order central gravitational attraction term. [Ref. 94:p. 26] 

The reentry vehicle is, however, considered to be in nearly vertical fall from 
an altitude of 30 km and below. Thus, there is only a minor dispersion of the impact 
point during the final seconds of flight. This rationale leads to the conclusion that it is 
not necessary to perform another transition from strongly perturbed Keplerian motion to 
an aerodynamic flight phase for the integration to Earth impact. Thus the Center Of 
Impact Window (COIW) is defined as the location at which the vehicle passes through 
the 30 km altitude. [Ref. 94:p. 26] 

The aerodynamic transition regime is defined as a computed, weighted mean 


Of Pusis ANd Pyss, Where 
P srs = f(h, >, A, UT, T, ly Bias Fio7 A,) (153) 


Puss = f(A) (154) 


which results in 
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(155) 


P = Wusis Pusis t (1 — Wass) Pusss 


where the weighting factor 
Wass (2) € [0,1] (156) 


is defined for the altitude region h,., = h = hg as 


Nona 
Wusis(A) = ze aes i 5 (157) 
where 
h = geodetic altitude 
6 = geodetic latitude 
rN = geographic longitude 
if = local solar time 


UT = universal time 
ta = day of the year 


When the altitude is high enough to maintain 


A 
Kn, = —>1 (158) 
d.. 
where 
d,. = characteristic length of vehicle 


OE, 


then Cp is equal to a constant given by C,“°, provided Kn, > 1. As the altitude 


decreases, a transition to hypersonic continuum flow is entered, where Kn, = f(h) and 


[Ref. 94:p. 27] 


Cy = Cy" + Ci" log, (C3 Kn.) (159) 
provided that 0.02 < Kn, < 1. During the hypersonic continuum phase another 
constant C, level is attained (about 50% less than the free molecular flow value), which 


is given by 


C,, = constant = Co” (160) 
provided 
Kn, = 0.02 
Ma, = 5 
The next phase of supersonic and transonic continuum flow can be 


approximated by an altitude dependent Mach number, Ma = f(h), where 


Cy = Cy” + Cy” (Ma, - 0.4) **C3” x exp|-C,"°(Ma, -0.4)] (161) 
provided 
Kn, = 0.02 
0.8 < Ma, <5 
Finally, in the subsonic phase, the drag coefficient is dependent upon viscous 


interactions and therefore Reynolds number, Re = f(h), given by 
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Coe Ci¥ (Re) **C," (162) 
provided 
Kn, < 0.02 
Ma, <= 0.8 
These dependencies, Cp = f(Kn, Re, Ma), with the underlying altitude 
functions Kn(h), Mach) and Re(h) are from the U.S. Standard Atmosphere (USSA) and 
are incorporated into the numerical reentry prediction software. These model constants 
are limited, however, to spherical and cylindrical shapes, with their longitudinal axis 
perpendicular to the airflow. [Ref. 94:p. 27] 
Since Ma = V,,,, / Vsoung and Re = V,,,, d/v are functions of geodetic latitude 
related by V,,,,a(h), kinematic viscosity, v(h), and aerodynamic velocity, V,,,,, the free 
fall of the reentering vehicle in the lower atmosphere may be determined by iteratively 


solving the following equation for the equilibrium descent velocity, V,,,, 


m_ gh) (163) 
A P(A) cp (A, Vin) 


Vz, = h? =2 
where 
g(h) = central gravitational acceleration of the Earth 
p(h) = local air density 


Using this technique, in the post-flight analysis of Salyut-7/Cosmos-1686, the 


following values are determined: 
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m/A = 159.5 kg/m? 


therefore V,,,, at 


h = 30 km, is equal to 330 m/s 
h = 20 km, is equal to 100 m/s 
h = 10 km, is equal to 70 m/s 


Figure 44 shows the ESOC reconstruction of the Salyut-7/Cosmos-1686 final 
descent altitude profile. The actual reentry occurred at 0345 UCT on 7 Feb. 1991 over 
(39.3°S, 69.7°W) [Ref. 94:p. 29]. Table 13 shows a comparison of reentry predictions 
from the U.S. and various European and former Soviet sites [Ref. 94:p. 32]. 

The conclusions of this European Space Agency (ESA)/ESOC report indicate 
that there was good correlation between ESA, U.S. and former Soviet predictions. There 
was some difficulty in interpreting the Soviet orbit determination results at approximately 
12 hours prior to reentry. Also at this time, the U.S. element sets became "time late,” 
due to transmission delay times from the U.S. to ESOC. Therefore, the ESA/ESOC 


predictions could not be maintained in "real time" after this point. [Ref. 94:p. 33] 


2. The LIFETIME Model 
The Aerospace Corporation previously developed the LIFETIME program 
and it was further refined through the work of reference [95]. This program is similar 
in some respects to the previously discussed ESA program, FOCUS. This program 
offers a "fast, efficient computer tool for orbital lifetime estimation” [Ref. 95:p. 17]. 
An advertised major benefit of the LIFETIME program is its usefulness as a “highly 


accurate, real-time reentry prediction tool." [Ref. 95:p. 20] This may also have 
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Figure 44: Salyut-7/Cosmos 


[Ref. 94] 
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Table Xill: SALYUT-7/COSMOS-1686 FINAL DECAY PREDICTIONS 
(Ref. 94] 












Re-Entry Predictions 
Prediction (times in UTC) 


Source Recelved Re-Entry Pred. 
at ESOC (COW) 


[NASA GSFC | O1-Feb 1458 | 06-Feb 21:99 
08-Feb 22:30 
06-Feb 04:00 
07-Feb 01:37 
07-Feb 03:38 
07-Feb 03:17 
CNES CST 06-Feb 18:23 
ESOC MAS 07-Feh 02:23 
07-Feb 03:05 
07-Feb 03:50 





CNES CST 


06-Feb 11:38 07-Feb 04:30 
CNUCE CNR 06-Feb 12:50 06-Feb 22:26 
06-Feb 14:09 07-Feb 03:28 


06-Feb 16:51 07-Feb 01:24 
06-Feb 16:58 07-Feb 04:14 


CNES CST 07-Feb 05:00 
06-Feb 04:19 
07-Feb 04:43 
CNES CST 07-Feb 04:38 
07-Feb 04:40 
07-Feb 04:05 


07-Feb 04:38 
07-Feb 04:20 
CNUCE CNR 07-Feb 04:00 

07-Feb 04:00 
07-Feb 03:57 
07-Fet 03:37 
07-Feb 04:29 
07-Feb 03.47 


07-Feb 13.00 07-Feb 03:45 


07-Feb 01:00 07-Feb 04:50 
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application as a first time modeling tool as compared to methods capable only of post 
flight analysis of reentry events. 

The LIFETIME program uses either the Jacchia-Walker (64) or Jacchia (71) 
atmospheric density models for altitudes above 90 km and uses the U.S. Standard 
Atmosphere (1962) for altitudes below 90 km. The latest version (4.0), the subject of 
reference [93], also allows differential corrections of the ballistic coefficient and 
movement of solar panels to simulate sun tracking. [Ref. 95:pp. 18-19] Version 4.0 was 
designed to solve a major deficiency of version 3.0, in that, there was a built-in 
uncertainty in Earth impact time prediction of at least one orbital revolution. This was 
inherent due to the minimum step size of one revolution, which was a function of using 
the averaged equations of motion. [Ref. 95:p. 22] 

A “unique” feature of the LIFETIME program is the differential correction 
of the ballistic coefficient. By using the least squares method, the following equation for 
differentially correcting the inverse ballistic coefficient (B*) can be formulated [Ref. 


95:pp. 28-29] 


(164) 





where 
N = number of observations 
a = a/a 
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= observed semi-major axis at time 1 
@ = semi-major axis at epoch 
e; | = observed eccentricity at time 1 

It can be shown that orbit decay in semi-major axis and eccentricity is 
directly proportional to the product of inverse ballistic coefficient and density (B*p). 
This differential correction process (which absorbs the atmospheric density uncertainties) 
results in a “converged ballistic coefficient" through multiple iterations of equation 175. 
This is how LIFETIME becomes more accurate than the atmosphere model it uses. [Ref. 
95:p. 29) 

LIFETIME integrates the averaged equations of motion using a step size of 
multiples of the orbital period, until the satellite reaches a specified or default altitude. 
Based on [86], the default altitude chosen is 42 nm (78 km) and this is where “breakup” 
is said to occur. LIFETIME then backs up to the last propagation step, resets the orbital 
elements and time to the previous step’s values. Next, the Runge-Kutta 7(8) integration 
routine is invoked and the satellite’s orbit is propagated through breakup to impact. (The 
elements are converted from classical mean to classical osculating elements and finally 
to Cartesian elements during the numerical integration process.) [Ref. 95:pp. 32-60] 
Numerical integration is the preferred method for dealing with the final stages of reentry 
to impact since it is the “most efficient and accurate way to handle regions where a high 


rate of change in the equations of motion" are prevalent [Ref. 95:p. 55]. 
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Several significant user input altitude values are: [Ref. 95:p. 59] 

1. RKALT - lower limit at which LIFETIME converts to Cartesian numerical 
integration. Default = 0.000 km, most effective if > BRKALT. 

2. BRKALT - vehicle breakup altitude. Default = 77.784 km. 

3. ENDALT - perigee altitude lower limit, ends integration and propagation. 
Default = 10.000 km, this should be set to 0.000 km to model impact. 

At the point where the satellite reaches the BRKALT, the lat/long projection 
of the satellite is recorded. This point is defined as the "debris heel point," or the first 
point along the groundtrack where debris could potentially impact. Propagation continues 
until ENDALT, where the lat/long is identified as the center of mass impact point. 
LIFETIME then resets the Cartesian elements to those at the BRKALT point, however, 
B is changed to 60 lb/ft? and the propagation to ENDALT is continued. This allows the 
computation of a "debris toe point," or the point of farthest travel of any breakup debris. 
The value of 60 lb/ft? for the ballistic coefficient is based on reference [86]. [Ref. 95:pp. 
59-60] Figure 45 shows the debris dispersion footprint as described above [Ref. 95:p. 
60]. Figure 46 shows the LIFETIME groundtrack and impact area plot as well as a 
sample altitude decay history [Ref. 95:p. 63]. (An item of interest not addressed by the 
author of reference [95] is the fact that the center of mass, in this sample plot, impacts 
after the "final debris impact.") 

The results of impact prediction accuracy comparisons between LIFETIME 
3.0 and 4.0 do not definitively show 4.0 as a significant improvement over 3.0. 


However, some notable improvement in accuracy for version 4.0 can be shown and 
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Figure 45: LIFETIME Debris Dispersion Footprint 


[Ref. 95] 
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Figure 46: LIFETIME Groundtrack And Altitude Decay History 
[Ref. 95] 
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greater confidence is held in version 4.0’s period calculations and final impact 
predictions. [Ref. 95:pp. 65-70] 


A sensitivity study and analysis was performed in order to determine the: 


1. Effects of NORAD data span on program accuracy 
2. Effects of prediction span on program accuracy 
3. Effects of solar conditions during differential corrections 


4. Effects of solar flux on overall program sensitivity 


The results of the sensitivity study and analysis are as follows: [Ref. 95:pp. 89-94] 


1. NORAD data span: 
(a) Period calculation accuracy - no strong trend for sensitivity noted 


(b) Impact prediction accuracy - general trend for less averaged impact time 
error for greater time spans 


2. Prediction span - impact time errors show a strong trend for increased accuracy 
with shorter prediction spans. 


3. Solar flux - the results indicated a complex relationship between the estimated 
solar flux for the prediction period, the actual values of solar flux, the 
prediction span, and the resulting impact errors. 

A general observation and conclusion is that fairly steady solar activity at 
modest levels provide for the most accurate impact prediction, especially when using the 


solar flux value of the last day on NORAD data as a constant during the prediction span. 


During periods of highly variable solar flux, the last day’s value of NORAD data may 
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not give the least error, however, this conservative assumption will not result in the 
maximum error either. [Ref. 95:p. 92] 
An appropriate concluding remark of reference [95] 1s __ that 


It seems there are still dynamic aspects of the final orbit decay and reentry process 
that warrant further study. [Ref. 95:p. 93] 


3. Modeling Ballistic Coefficient 
Analysis of historical data shows that in the final days prior to decay, the 
ballistic coefficient, B, (CpA/m) exhibits a pronounced variation with time. This 


variation is due to one or more of several causes: [Ref. 92:p. 2] 


1. Tumbling or weather vaning of the satellite 
2. Loss of mass 


3. Variation of the drag coefficient 


Presently, none of these effects are modeled explicitly since this requires precise 
knowledge of the multiple variables described in Chapters II and II, such as vehicle 
mass, attitude and shape. Additionally, any unmodeled forces acting in the in-track 
direction, such as those inaccuracies in the atmospheric density model, can manifest 
themselves in the solution ballistic coefficient. [Ref. 92:p. 2] 

Despite the lack of specific knowledge of unknown forces acting on a 
decaying satellite, an improvement in reentry prediction accuracy can be made by 
introducing a new model parameter. This parameter takes into account the fact that for 


many decaying satellites, the variation in B appears to be linear during the last several 
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days prior to reentry. This new parameter is By-dot, the time derivative of the ballistic 


coefficient, which 1s parameterized as 


B = B, + B,(t-t) (165) 
where subscript 0 indicates the quantities at the solution state epoch, ty. [Ref. 92:p. 3] 
The introduction of By-dot requires that new semi-analytic partial derivatives 
for By and B,-dot be calculated. These must be determined over a finite time interval 
using a differential correction fitting process. Several such differential corrections must 
be made over the course of the final days of orbit decay and reentry, with a new By and 
B,-dot being determined for each differential correction. It is necessary to start the 
fitting process with close approximations of B, and B,-dot for the first differential 
correction. These initial "guesses" come from analysis of historical data. [Ref. 92:p. 3] 
Reference [92] studied 264 reentry events over the time span from Jan. 1977 
to Mar. 1979. All of these objects were processed as TIP reentries and historical records 
from the SSC were used in this investigation. In order to optimize the data, the records 
of 156 reentries were chosen, since these represented rocket bodies and payloads from 
the former USSR, which presumably would have well-known physical attributes. Table 
14 shows the data distribution used throughout this investigation. [Ref. 92:pp. 3-4] 
The development of the solutions for the partial derivatives of By are not 
shown here for the sake of brevity, however, the results are as follows [Ref. 92:pp. 10- 
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Table XIV: DISTRIBUTION OF DECAYED ORBITAL OBJECTS 
(Ref. 92] 


DISTRIBUTION OF 264 TIP LOG DECAYED ORBITAL OBJECTS 


YEAR: 1970 

MONTH: JF MA MJ J A S$ ON D 
NUMBER: 2 

~ YEAR: “1997 - 
MONTH: J F MAM JI 3 A S ON D 
NUMBER 8 8 116 158 8 11149 8 7 
YEAR: 1978 

MONTIL: JF MAM JI J A S ON 
NUMBER: 119 9 13 6 127 8 11:13 8 18 
YEAR: 1979 

, MONTH: J oF MA J J A S ON D 
NUMBER; 6 10 7 

YEAR: 1983 

MONTH: Jo F.M A MJ 43 A S ON D 
NUMBER: 1 


Swe es Gee eeSOs eG OG GOe S@2 Fee Ge Gee eOesTaevsteeweoGOnaeeest 2G eGGasteeeensewvwoeeoeseecoeneaaeanenezea es 2 ®@ 


© Selected 156 Soviet orbitally decayed object distribution:§.7.8 


o ROCKET BODIES DIA(m) LNG (m) 

21 A2 VOSTOK ( 440 kg) 2.6 7.5 
89 A2 VOSTOK (2500 kg) 2.0 2.0 

4 Bl Cosmos (1500 kg) 1.65 8.0 

8 Cl SKEAN (2000 kg) 2.0 6.0 

3 Dl PROTON (1900 kg) 3.9 4.0 
3  #Di- PROTON (4000 kg) 4.0. 12.0 
-3 FIM SCARP (700. kg) - 2.0-: 5.0 

1 FIM SCARP (1500 kg) .:, 2.5 6.0 
PAYLOADS$ . 

3. NAVSAT ( 875 _ kg) ? ? 

1 ( 325 . kg) 

1. (3640 kg) 

| ene ( 680 , kg) 

8 MILSAT. ( 400” kg) 

Z MIL/SCI ("550° kg) 

4  RORSAT® (4500 Kg)° 

2 MIL ( 900 - kg): 

1 (1120 kg) 

T (3800,. kg). 


§° The exact. size and shape of the payloads are unknown at this point. 
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da 1 
eae See 
B, a By ayy ie 

da 
B,— = B, 60 rle 
6, % Be 


= -—B, —]| —adt 
"0B, 2 °aln B 
OL 3 0 
— =- --—B,—(t-t —adt 
°0B, 4 a oJ 
p&p oF 9 
OB, OB, 
where 
ar =ecos (2 + w) 
a  =esin (22 + w) 
é = time rate of change of eccentricity due to drag perturbations 
a = time rate of change of semi-major axis due to drag perturbations 
L = mean longitude 
X = tan (1/2) sin 2 
Y = tan (i/2) cos 2 
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(166) 


(167) 


(168) 


(169) 


(170) 


Similarly, the partial derivatives of By-dot are given [Ref. 92:p. 12] 


B, = = By ay flew 5 ear (171) 

= = B, a0 J let) 7 =a (172) 

8 = Boag f|(t-W 5 ala (173) 

or - “Gan (t-t,) fije-w 5 ala (174) 

Bok = BY = 0 (175) 
aB, dB, 


Of the 156 cases studied in the development of the new parameter and its 
associated partial derivatives, two case studies were chosen for the application of the new 
method. These cases were Cosmos-954 and Cosmos-1402, both of them being significant 
due to their potentially survivable nuclear reactors. [Ref. 92:p. 13] 

The conclusion of this paper is that the new method is more accurate than the 
conventional method as evidenced by the comparison of the best fit differential correction 
root mean square (RMS) of residuals for each of the 25 sliding-fit runs over the last nine 


days of orbit for Cosmos-954. The same comparison 1s made for Cosmos-1402 with the 
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only exception that 26 vice 25 SP differential corrections were run. The results of these 
comparisons are shown in Figure 47 [Ref. 92:pp. 15, 17]. 

As Stated previously, the historical data analysis of most satellites showed the 
linear trend in B only over the last 3-5 days prior to reentry. However, one of the two 
specific cases chosen for application of the new method, Cosmos-954, did not exhibit this 
trend until the final day prior to reentry. Table 15 shows a comparison of decay time 
error root mean square over the last five days of orbit (for these two case studies) as well 
as a similar comparison after the linear trend in B becomes apparent in both reentries. 
[Ref. 92:p. 18] 

Another method using a correction of the ballistic coefficient comes from 
reference [47]. This investigation states that during the joint international effort to 
predict the reentry/impact of Salyut-7/Cosmos- 1686, it was observed that the "most exact 
data" had consistently been presented by the Space Control System (SCS) of the former 
USSR. This was explained by the fact that: 

...great attention has been given to the problem of determination and the prediction 
of the satellite motion in the atmosphere: for the solution of this problem they have 
been conducting the cycle of theoretical and experimental studies during [the last] 
20 years. [Ref. 47:p. 35] 

The model used by the SCS to describe satellite motion will, to a great 
extent, determine the "completeness" and precision of the time and area of low-orbit 


satellite reentry predictions. Gravitational field models in spherical functions of 36 or 


higher orders are not necessary, they only serve to increase the computation time. At 
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Figure 47: Best Fit DC RMS For Cosmos-954/1402 
[Ref. 92] 
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Table XV: DECAY TIME ERROR RMS 
[Ref. 92] 


—=—=—=— —— SS 7 


DECAY TIME ERROR ROOT MEAN SQUARE (RMS) OVER LAST 5 DAYS 


Decay Error RMS(hrs) 
Decay Error RMS(hrs) Bo, By Solution, 


Satellite Bo Solution By=0.0 for Predictlon [mprovement(%) 
COSMOS 954 3.51 | 4.11 -17.09 
COSMOS 1402 Nees; 1.43 10.06 


DECAY TIME ERROR ROOT MEAN SQUARE (RNS) AFTER B BECOMES LINtzAi¢ 


Decay Error RIMS(irs) 
Decay Error RMS(hrs) Bo, B, Solution, 





Satellite Bo Solution B =0.0 for Prediction Improvement(%) 
COSMOS 354 1.22 0.48 60.66 
COSMOS 1402 1.51 0.98 35.10 
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low altitudes, two to three days prior to reentry, a simple geopotential model is suitable, 
since it is the atmospheric resistance force which is increasing by several orders of 
magnitude, not gravitational forces. [Ref. 47:p. 36] 

Prediction of low-orbit satellite motion is accomplished by numerical and 
semi-analytical integration methods. The semi-analytical methods are based on the 
asymptotic solutions given by Knilov-Bogoljubov, Zeipel-Brouwer averaging methods. 
[References given are written in Russian and the authors were unable to obtain English 
translations.] [Ref. 47:p. 36] 

The actual method applied to a reentry event may be a variant of some 
numerical and semi-analytical methods. These specialized fast-acting, semi-analytical 
prediction algorithms for low-altitude orbits have “counting times of order 0.01 seconds" 
for a 24 hour prediction span. [Ref. 47:p. 36] 

The reentry of the Salyut-7 complex was predicted and confirmed by 
observation in the range of 75-105 km. In this altitude range, it was also observed and 
predicted that breakup would occur. [Ref. 47:p. 39] 

It was assumed that the ballistic coefficient would vary during reentry, to the 
same degree as it did prior to an attitude control maneuver. This assumption allowed the 
reentry problem to be solved by varying the initial conditions within their ranges of 
possible errors, with subsequent integration of motion equations up to the point of 
reentry. In many cases it was only the ballistic coefficient which must be varied. This 
sensitivity of reentry time to the ballistic coefficient is given by the following 


approximate formula [Ref. 47:p. 39] 
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Vie a oe (176) 


where 
t; = time interval of motion prediction to the moment of reentry 
k = ballistic coefficient = CpA/2m 


Given a 1% change in ballistic coefficient, the down range travel is altered 
by 180 km. A 2-3% change in ballistic coefficient (as seen prior to the attempted control 
maneuver) alters the predicted impact point by + 500 km. Figure 48 shows the 
dependence of time, latitude and longitude on the ballistic coefficient [Ref. 47:p. 40]. 

The approximate estimate of the probable debris fallout range, considering 
the errors in determining the reentry point, with the vehicle destruction occurring at 75 
km altitude, is equal to + 1000 km relative to the calculated impact point. This 


corresponds to the sub-satellite track crossing Chili and Argentina from the south-west 


to north-east. [Ref. 47:p. 40] 
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Figure 48: Dependence Of Time, Latitude And Longitude On B 
[Ref. 47] 
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C. FACTORS INFLUENCING REENTRY/IMPACT DISPERSION 
The reentry vehicle’s deceleration may be expressed in terms of "gravitational 
acceleration forces” (g’s) as follows [Ref. 96:p. 1] 


1 p y? 
Dynamic pressure _ 2 


Ballistic coefficient W y 


(177) 


Cy 


Assuming that density and wind velocity are the two atmospheric characteristics 
which most strongly influence the motion of a reentry vehicle passing through it, 


reference [96] states the following: 


1. Reentry vehicle is non-ablative and non-maneuvering (W and A are constant) 
2. Drag coefficient is dependent upon altitude (Mach dependence) 


3. Drag coefficient is shape dependent 


The displacement in impact point can be considered to be the sum of miss 
contributions from the various atmospheric layers, given by the product of influence 
coefficients times changes in density or wind velocity for the respective altitude layers. 
[Ref. 96:p. 2] 

Reference [97] studies the effects of a non-spherically symmetric atmosphere on 
the reentry/impact location of a decaying satellite orbit. This investigation deals with two 
specific cases, a 90 degree and 30 degree inclined, low Earth orbit and shows the effect 


of the atmosphere’s diurnal bulge on the impact location. It is observed that a decaying 
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object is "more likely" to reenter and impact the Earth at certain latitudes than at others. 
[Ref. 97:p. 111] 

The principal effects contributing to this impact distribution, in addition to the 
diurnal bulge, are the geometric oblateness of the atmosphere and the gravitational 
oblateness of the Earth. Also identified as important parameters, are orbital inclination 
and ballistic coefficient. 


The results of this investigation are as follows: [Ref. 97:p. 43] 


1. Polar orbit (90° inclination) reentry is most likely in the equatorial latitudes. 


2. Latitudes of maximum or minimum likelihood for impact are not significantly 
affected by changes in ballistic coefficient. 


3. Large, "balloon type” (spherical) satellites react to the diurnal bulge by causing 
impact on the opposite side of the Earth. 


4. Smaller, heavier satellites show little influence of the diurnal bulge, but they 
exhibit a greater variation in impact probability with latitude. 


5. The impact distribution for a given spherical satellite "flattens" as the inclination 
is moved from 90° to 30°. This is attributed to the fact that the extremes in 
Earth radius differ by only 5 km under the 30 ° inclined orbit as opposed to 21 
km under the polar orbit. 


Reference [98], a survey of satellite lifetime and orbit decay prediction, conducted 
in 1980, states: 
...gravity perturbations cannot, by themselves, lead to orbital decay since they are 
conservative. But, they can induce oscillations in the orbit...drag is proportional 


to atmospheric density...gravity perturbations coupled with drag are more 
significant than when only drag is modeled. [Ref. 98:p. 3] 
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Also noted is the observation that non-symmetric mass distribution and shape may 


amplify the torques created by: [Ref. 98:p. 4] 


1. Gravity gradients 
2. Aerodynamic forces 


3. Solar pressure forces 


An interesting note, also, is that any mass loss increases the magnitude of the drag 


deceleration since [Ref. 98:p. 5] 


Bag — (178) 


A final noteworthy observation from reference [98] is that any attempt to model 
satellite attitude, shape, lift, mass distribution, ablation and breakup would largely exceed 
the level of sophistication currently available (1980) for atmospheric density models, 
therefore, these non-atmospheric factors are neglected and treated as higher-order effects. 
[Ref. 98:p. 5] 

Reference [99] was initiated after the reentry of part of Sputnik IV over Wisconsin 
in September 1962. [Ref. 99:p. 2] This study of factors influencing the prediction of 
orbit decay and impact states: 


...final reentry occurs shortly after the point of minimum altitude (not coincident 
with perigee, because of the Earth’s oblateness) [Ref. 99:p. 29] 


It is possible, knowing only the initial inclination angle and reentry vehicle orbital 


path, to determine the points at which the minimum altitude points will occur. 
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Therefore, knowing where these points occur, it may be possible to predict reentry within 
several tenths of an orbit revolution, as opposed to ambiguities of one complete 
revolution as noted in reference [95] (LIFETIME v3.0), since reentry tends to follow the 
point of minimum altitude. [Ref. 99:p. 33] 

Without such precise knowledge of the body’s ballistic coefficient, the 
characteristics of the atmosphere, the orbital elements and their variation, only 
"probability estimates” can be made of impact latitude, based simply on the relative time 
spent by the satellite in various latitude bands. This conclusion yields the following 
relationships: [Ref. 99:p. 33] 

1. If the initial inclination equals 50° - three times more likely to impact between 
40-50° band than 0-10° latitude band 

2. If the initial inclination equals 49° - two times more likely to impact the 
continental U.S. than a 65° inclined orbit 

A final comment about factors influencing reentry and impact prediction from this 
investigation is that: 

In spite of the present reentry rate of about one object a day (most of them burning 
up, and 2 out of 3 of Russian origin), the threat imposed on the Earth population 
from direct hits by debris from uncontrolled reentries of space objects can 
generally be regarded as minor when compared with other risks of civilization 
(traffic accidents, etc.). An accumulation of worst case assumptions could lead to 
one casualty every five years for a densely populated area of the size of the Federal 
Republic of Germany in case of no prior warning. In case an early warning could 
be issued, such casualties could most likely be avoided completely. [Ref. 100:p. 
49] 


Figure 49 shows the magnitude of the low Earth orbit problem purely as a function 


of the number of objects in orbit [Ref. 100:p. 51]. Figure 50 shows the limiting factor 
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Figure 49: Near-Earth View Of All Cataloged Space Objects (1987) 
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of geographic location of observation sites from which data may be obtained for 
reentry/impact predictions [Ref. 100:pp. 5-6]. 

The final factor in reentry/impact prediction to be discussed here is the rotation of 
the orbital plane during the final phase of reentry, as presented in reference [101]. This 
study states that planar motion is representative of, or an adequate model for, reentry 
deceleration and heating rate studies. Itis, however, inadequate for the study of ballistic 
reentry below altitudes of 60 km and ballistic coefficients in the ranges 0.001 m’/kg < 
B < 0.1 m’/kg. For these values of B, the orbital plane begins to rotate at altitudes of 
30 km and 60 km respectively. [Ref. 101:p. 1215] This investigation deals exclusively 
with the final descent from 120 km, and defines this as the "reentry phase” [Ref. 101:p. 
1216]. 

As previously mentioned, maximum heating rates and deceleration occur at altitudes 
high enough, such that, rotation of the orbital plane is negligible in the calculation of 
these parameters. However, in the calculation of the total range traveled to Earth 
impact, these lower altitude effects and orbital plane rotations become significant. [Ref. 
101:p. 1217] 


Several initial assumptions are made: [Ref. 101:p. 1217] 


1. 0.001 m’/kg <= C,A/m < 0.1 m’/kg 
2. O < Wy <= 45° , where Wy is the initial reentry trajectory angle 


3. Vo <= 15 km/sec 
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As the vehicle descends, the ambient air becomes denser and the rotating 
atmosphere begins to force the vehicle’s motion towards a direction parallel to that of the 
atmosphere’s velocity. This effectively rotates the vehicle’s orbital plane towards the 
direction of the atmospheric velocity vector. Figure 51 shows this rotational effect on 


the orbital plane [Ref. 101:p. 1219]. 


z-axis points to the north, r-axis to the east, the plane + is the surface of the ‘inertial Earth’ 





Figure 51: Orbital Plane Rotation Due To A Rotating Atmosphere 
(Ref. 101} 
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V. STOCHASTIC AND STATISTICAL PREDICTION METHODS 


This chapter will describe three investigations that used stochastic and statistical 
techniques and methods to predict the time and impact location of decay-induced 
reentry. Two investigations have direct applicability to the current reentry model used 
by the USSPACECOM and the third investigation is a Monte Carlo simulation analysis 


that was used to predict the footprint dispersion area of Skylab. 


A. EXTENSIONS OF THE PHYSICAL MODELING REENTRY THEORY 

In Chapters II and IV, the current reentry theory used by the USSPACECOM to 
predict impact location and time was presented. The theory’s accuracy is limited by the 
ability to observe the reentry process and the inherent deficiencies in the model used to 
represent the physical reentry process. This section will examine a Doctoral dissertation 
and a Master’s thesis from the Air Force Institute of Technology (AFIT) which are 
proposed extensions or modifications of the current model. [Refs. 48,102] These 
investigations utilize stochastic and statistical methods and techniques to improve the 


predicted impact location and time. 


1. Estimation of Reentry Trajectories 
Reference [102] developed a linearized differential corrector technique as an 
extension of the existing orbital estimation technique into the reentry region to estimate 


the reentry trajectory. Reentry observations from a space based sensor capable of 
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providing infrared angular trajectory observations at fixed, discrete time intervals and 
large uncertainties in true reentry dynamics, were the primary engineering considerations 
of this problem [Ref. 102:p. 1]. The author states the following about the current 
deterministic model: 
...the basic differential corrector with deterministic dynamics is inadequate for 
processing true reentry data. The deterministic dynamics and infinite memory of 
the basic formulation, cause the estimator to yield significantly biased solutions 
relative to the standard deviations of the estimator-computed covariance matrix. 
These occur when processing data reflecting the dynamic variations anticipated in 
true reentry trajectories. The use of the estimator in this form would not provide 
an accurate estimate of Earth impacts for satellite debris. [Ref. 102:p. 19] 
Since the uncertain dynamics of the reentry process pose many difficulties for the 
existing deterministic model, a significant portion of the research in the dissertation was 
devoted to identifying limits of various estimation techniques that were considered for 
this problem. [Ref. 102:p. 1] 

An adaptively determined, ad hoc, scalar finite memory or fading memory 
parameter approach was selected. The motivation of this approach is based on the fact 
that the existing reentry theory attempts to group all unmodeled parameters as a constant 
with the ballistic coefficient (C,>A/m) over some short trajectory span. [Ref. 102:p. 11] 

In order to extend the existing orbit determination technique into the satellite 
reentry problem, reference [102] defined a reentry dynamics model and developed the 
weighted least-squares differential corrector structure currently used in the deterministic 


model. The estimator update and propagation equations of the differential corrector 


structure for an infinite memory formulation are summarized as follows: 
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function of the observation residuals. 


Propagations between epochs: 
State: integrate with initial conditions x,(t,,), from the epoch at t,, 


x(t) = fix(D,) 


to obtain x,(t, 4 ,) at epoch t, 4, _ 
Covariance: 


Scat OUaeish) Sone Ceeneeem 


date at nex h: 
State: 


I 
5 +1, + J zz X(t + y Bi » x(t, + Vj 


i=] 


Covariance: 
_ =I a alga 6 
Aree Ome Sy a Ley 


(179) 


(180) 


(181) 


(182) 


Recall that 1 is the number of iterations required to satisfy convergence and T,) * 
Ray Tg is evaluated from the reference trajectory on the final iteration, 1. [Ref. 


102:pp. 33-34] 
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Equation (179) is a nonlinear dynamics model, where f(x(t), t) is a deterministic function 
of the state variables and is a continuous function of time. The variable, 6x(t), from 
equation (181) is a "seeked" correction that will minimize a weighted quadratic cost 
Iterations of the differential corrector continues 
until the observation residuals convergence criteria is satisfied. Equation (180) is the 
initial state covariance matrix that remains constant until the iterative process has 
converged. The quantity $(t,,,,t,,) from equation (180), is the state transition matrix. 
The final iteration in the process yields equations (181) and (182). Rand T,) from 
equation (182) are respectively the matrices of the observation noise covariances and the 


observation index being processed. [Ref. 102:pp. 20-29] As previously mentioned: 


Application of this batch processing algorithm to reentry estimation has often 
resulted in poor estimator performance. This is largely due to the more significant 
non-linear dynamics of reentry and the use of a deterministic and simplistic 
dynamics model ... in an uncertain dynamics region. [Ref. 102:p. 31] 


In the development of the reentry dynamics model, the author chose an eight 


dimensional state vector defined as [Ref. 102:pp. 41-42] 


ei uk 
x, =% 
x3 a 
x,=y 

x(t) =| (183) 
XK, = 2 
X, = 2 
x, = Bpy 
x, = Q 

where 
%,y¥,Z = velocity components (184) 


X,Y,Z, = position components (cartesian) 


B = ballistic coefficient (CpA/m) 

Po = Sea level air density from equation (10) 

Q = density scale height (RT) / Mj& ) from equation (10) 
M, = molecular weight of air 

Zo = acceleration of gravity at Earth surface 

T) |= atmospheric temperature at sea level 

R = gas constant 
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H’ = altitude (g/g ,H,) from equation (10) 


gZ = local gravitational acceleration from the geopotential model x,y,z 
coordinates 
2) = reference geoid level gravitational acceleration 


H, = geocentric altitude (R-R,) 
R = local radius position relative to the Earth center (x? + y? + z’)!” 


R, =R(1-f[l - 2f- f) cos? 6] !” 


fe = flattening factor of the reference geoid whose elliptical shape is consistent 
with J, 

5 = local latitude ( cos'[ (x? + y?)!” / Ro] ) 

Ry = radius of the reference geoid at the equator 


The exponential atmospheric density model, equation (10), was used in the 


reentry estimator because of: [Ref. 102:p. 37] 


1. The reduced mathematical complexity. 
2. The availability of continuous valued density and partial derivatives of the 
density along the reentry trajectory. 

The eighth state variable, Q, was chosen since it is slowly changing at the 
reentry altitudes ( < 100 km). For this application, the initial value, Q = 7.0031 km 
and its covariance were derived using a least squares fit to the base density values of the 
altitude layers from the U.S. Standard Atmosphere 1962. Also, an Earth Centered 
Inertial (ECI) coordinate frame was chosen to minimize the complexities of the 


observation noise covariance matrix, Rq. [Ref. 102:p.40] 
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Both aerodynamic and geopotential terms within the deterministic function of 
State variables, f(x(t) , t), are included in the estimator dynamics model. The dynamics 


equation is given in the form [Ref. 102:pp. 42-44] 


a, = % 

4, = Sa, * Sp, 

eX, 

a2 x, = Sa + Sp (185) 

os = *6 

3 0 eae 

X, = 

i, = 0 


The aerodynamic acceleration x,y,z components from equation (185) are 


derived from equation (6) 


iis -=Bp V, &+ oy) (186) 
fa, = ~5BeV, (0-0) (187) 
f, > “5B V,2 (188) 


The geopotential acceleration terms from equation (185) were derived from 
the Smithsonian Astrophysical Observatory SAO-III Earth Model using the model 


parameters and zonal and tesseral coefficients 
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f, = cosa,VV, - sina,VV, (189) 


fi, = sina VV, + cosa, VV, (190) 
i = VV. (191) 
where 


a, = Right Ascension (RA) of the Greenwich Meridian 


The gradient terms from equations (189) through (191) are defined as 








=e aU.” ov.” 
VW, = > > (c,, dae are Sa 2 | (192) 
n=O m=0 








= 4 C ay > (193) 








= ee ou,” ov,” 
mB EM os 


where 
Cin = zonal harmonic coefficient 
Sin = tesseral harmonic coefficient 
M R; 


Uy = 2 Pr) cos(ma) (195) 
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_ GM Ry 





ye —y P,"(3') sin(ma) (196) 
G = universal gravitational constant 

M = mass of Earth 

Ry |= mean equatorial Earth radius 

R = distance from Earth center 

P,” = associated Legendre polynomial of degree n and order m 

6° =2/R (z = coordinate) 

r = longitude 


As previously mentioned in Chapter II, in the region below 120 km, the first-order zonal 
harmonic, J,, may attain a magnitude approaching that of the atmospheric drag. Based 
on this fact, in the development of the actual estimator and simulation, only the central 
gravity, Co, and Earth oblateness (J,), Cy) , terms were used in the reentry altitude 
regions. 
Numerical simulations using simulated reentry data derived from a realistic 

"truth model" were conducted on the basic estimator structure and dynamics model to 
quantify its performance. 

These analyses examined the effects of mismatch between the truth model and the 

estimator model dynamics, accuracy variations on the angular observations, 


multiple orbital observation locations and variations of the geometric relationships 
between the observing satellites and the reentry trajectory. [Ref. 102:p. 14] 
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A series of Monte Carlo analyses were conducted which showed the model dynamics as 


the most important factor impacting the estimator performance. This was accomplished 


by: 


1. Considering the velocity estimates and mean bias of the trajectory position vs 
the magnitude of the standard deviations from the estimator computed state 
covariance matrix. 

2 Comparing the magnitudes of the standard deviations from the estimator 
computed state covariance matrix to those derived in the Monte Carlo samples. 

In the search to find possible solutions to the problems associated with the 
estimator performance, reference [102] reviewed and discussed the limitations of several 


model compensation methods that are relative to the reentry application, including: [Ref. 


102:p. 15] 


1. Adding a pseudo-noise compensation to the model dynamics 
2. Adaptive estimation methods 


3. State covariance deweighting techniques 


The author developed a fading memory differential model compensation method based 
on the fact: 


...the previous numerical results indicated, the fundamental limitation of the infinite 
memory estimator formulation with deterministic dynamics is the biased estimator 
solutions which occur when processing true reentry dynamics. With 1) exact 
dynamics 1i) an upper limit on the time span of valid linearization, and 111) a 
lower limit on observation accuracy (greater than or equal to 10° radians), 
acceptable estimator performance is available in terms of bias and RSS/(ONE 
SIGMA) Ratio. [Ref. 102:p. 112] 
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Since the variations in global atmospheric density and the dynamic changes 
of the reentry body pose such a virtually intractable problem within a deterministic set, 
improving the linearized estimator performance is based on the impact of uncertainties 
within the ballistic coefficient and the atmospheric density [Ref. 102:pp. 112-113]. 
Improvement of the estimator was achieved by an adaptive determination of an ad hoc 
scalar multiplier, y. A finite memory on the processing of earlier observations 1s 
implemented by multiplying the parameter, +, to the terms of the state covariance matrix 
prior to an observation update. At each update along the trajectory, the size of the 
change in state variable, 6x; 1s examined. Each 6x; magnitude is then compared to the 
magnitude of the standard deviation of their respective terms within a "deweighted" 
covariance matrix which 1s computed by the estimation algorithm. As with the infinite 
memory estimator model, a series of Monte Carlo analyses were conducted to assess the 
model’s ability to estimate an anticipated true reentry dynamics trajectory. [Ref. 102:pp. 
15-16] 

Finally, reference [102] demonstrated the ability of the fading memory 
method to provide a tangent plane projection of Earth impact locations as shown in 
Figure 52 [Ref. 102:p. 192], by using bias magnitudes within the standard deviations of 
the deweighted state covariance matrix. 

The standard deviation of the position error from the deweighted state covariance 
matrix provides a good definition of the uncertainties in the estimated Earth impact 
locations, thus it can be used to define a search area for recovery of satellite 
debris. These results illustrate the viability of the method to estimate decayed 


satellite impact locations and uncertainties significantly improved over existing 
astrodynamic applications. [Ref. 102:p. 16] 
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Figure 52: Tangent Plane Projection Of Earth Impact Location 
[Ref. 102] 
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The primary conclusions and recommendations from the dissertation are: 


[Ref. 102:pp. 198-202] 


1. Dynamics uncertainties of the general satellite decay trajectories significantly 
affect the estimator performance as shown by the Monte Carlo simulation 
analyses. 


2. Angular observation accuracies with standard deviations less than 10° radians 
induce significant error in the state estimate vs the standard deviations of the 
state covariance matrix in deterministic dynamics models. 


3. A recursive formulation of the estimator is recommended that uses a short time 
span between the epoch or trajectory update point and the observation(s) being 
processed due to the anticipated dynamics uncertainties of reentry. 


5. Multiple observations from more than one the orbital source are required to 
improve the observability of the reentry and to provide higher data content over 
similar time spans in-order to achieve acceptable estimator performance in terms 
of bias and RSS/(ONE SIGMA) ratio. 


6. An eight dimensional state vector provided superior estimation performance as 
compared to the seven dimensional vector used in the current deterministic 
models. Performance improvement is achieved through simpler mathematics in 
the dynamics model and continuous valued partial derivatives of the dynamics 
over the trajectory space for a Taylor’s series linearization. 


7. The adaptive,ad hoc scalar fading memory parameter is easily incorporated into 
the basic estimator structure. 


8. A Monte Carlo derived impact covariance for the final propagation phase to 
impact preserves the integrity solution statistics over the non-observable final 
portion of the trajectory. Impact location uncertainties are on the order of one 
to two magnitudes smaller than those available from current operational 
techniques. 


9. Further investigations which vary the observational data rate and the time 
varying character of true reentry dynamics are needed to examine the estimator 
performance extensions. Specifically, more accurate observations (much less 
than 10° radians) with higher data rates and frequency variations as well as 
alternative measurements, such as range or range rate. 
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10. Further analysis of applying the fading memory method to very high altitudes 
and shallow reentry angles is recommended due to the fact that violent dynamic 
changes under these conditions may result in divergent estimator performance. 

11. Application of this estimation technique to a wide class of reentry trajectories 
could provide a large empirical data base to improve the estimator dynamic 


model. Dynamic model pseudo-noise compensation, statistical linearization, or 
higher order filters investigations should be pursed. 


2. Analysis of Tracking and Impact Prediction (TIP) 
Reference [48] analyzed the accuracy of early TIP processing (Chapter II, 
p.58) conducted by the USSPACECOM on 180 objects that decayed during the years 
1987 through 1990. As part of the analysis, early TIPs (7 day to 3 hour predictions) 
were compared to the final TIP. The time error for each TIP run was calculated and 
compared to the + 20% accuracy level claimed by the SSC at Cheyenne Mountain AFB. 
Results from the comparison study indicate: 
1. The decay prediction accuracy is usually, but not always within the claimed 
accuracy level as shown in Figure 53 [Ref. 48:p. 37]. 
2. The existence of a positive bias which indicates early TIPs are routinely late 
relative to the final TIP. 
The author developed six multiple linear regression models that could be incorporated 
into the TIP decay procedures in an attempt to model out some of the positive bias found 


in TIP decay prediction data. [Ref. 48:p. viii] 
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Figure 53: Decay Predicted Accuracy (By Year) 
[Ref. 48] 
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The data used in the investigation was collected from: [Ref. 48:pp. 15-17] 

1. TIP Required Item Checklist--a manually kept chronological account of each 
TIP processing. 

2. Decay History--a computer generated log of each SP update that includes the 
run time, time of last observation, epoch time, epoch revolution number, B-term 
(ballistic coefficient), period and decay prediction. 

3. Final Element Set--a listing of the final orbit parameters. The eccentricity and 
mean motion data from the listing were used in the investigation. 

Only those objects that received the entire TIP update cycle (7 day through 
final run) were used to accurately analyze the effects of each successive update and 
prediction [Ref. 48:p. 17]. In order to assess the accuracy level claimed by the SSC, 
reference [48] used the data from the final decay prediction as the control by which to 
compare the earlier predictions. 

The final prediction was chosen as the control because it uses observations which 
are closest to the actual impact point and is considered to be the most accurate 
prediction available. In order to further justify the use of the final prediction as 
the control, a statistical analysis was also performed to directly compare the few 
sighted reentry points, called Vis Obs, with the final prediction made by the Space 
Surveillance Center. [Ref. 48:p. 18] 

Of the 180 objects studied, 93 were Vis Obs. Figure 54 shows the mean time 
error of the final run time vs the Vis Obs [Ref. 48:p.25]. The size of the final time error 
standard deviation decreases as shown in Figure 55 [Ref. 48:p. 26]. This decrease may 
be correlated to the level of solar activity during the 1987-1990 time period. Solar 


activity levels began to dramatically increase in 1987 and continued to increase through 


the solar maximum (March 1990). The rate of change increase of the sunspot activity 
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Figure 54: Final Run vs Vis Obs Mean Time Error 
[Ref. 48] 
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Figure 55: Final Time Error Standard Deviation 
[Ref. 48] 
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and solar flux during 1987-1988 was greater than in period from 1989-1990. According 
to the author, the trend shown in Figure 55 may be related to the lesser rate of change 
during 1989-1990. [Ref. 48:pp. 26-28] 

The time error for each separate TIP run was calculated as the difference 
between the predicted decay time for that run and the final run [Ref. 48:p. 18]. Figure 
56 shows the graphic calculation results for the mean time error and time error standard 
deviation [Ref. 48:p. 29]. 

In addition to the time error calculations, the location error was calculated 
by taking the difference between the predicted location point for each run and the final 
predicted point. The method used the mean motion (n) from the Final Element Set to 
accurately determine the velocity for each TIP object. Since mean motion was 
unavailable for the early TIP runs, an approximation was made using the final mean 
motion value for each TIP run. This introduced an error into the calculation of the 
location error. However, because the location errors are very large (thousands of km), 
the error introduced by using the final mean motion instead of the actual mean motion 
for that particular TIP run was considered insignificant. The mean motion value was 


used to calculate the semi-major axis distance from the equation 


1 
7 eal (197 


The velocity was calculated using the semi-major axis from the equation 
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Figure 56: Mean Time Error And Standard Deviation (1987-1990) 
[Ref. 48] 
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— »|2 : 4) (198) 
r a 

The object’s velocity was then multiplied by the previously calculated time error in order 
to determine the location error (difference between the predicted location point for each 
run and the final predicted point). [Ref. 48:pp. 19-22] Figure 57 shows the graphic 
calculation results for the mean location error and location error standard deviation [Ref. 
48:pp. 32-33]. 

Multiple linear regression describes the relationship between several 
independent variables and a dependent variable. The motivation for developing the 
multiple linear regression model to eliminate the bias found in the mean time error was 
based on the second objective of the investigation. This objective was to determine if 
it was advantageous to initiate an OPREP-3 report (used to notify higher authority of a 


potential reentry within 100 nm of the former Soviet Union border) earlier than the 6 


hour mark. The first step in the process was to determine if the model in the form 


E(t) = By + Bit, + Bot, + Baty + Buty + Bots + Bots + Bat (199) 
where 
E(t,) = expected value of the final decay prediction time 
t = early TIP decay predictions 
te = final decay prediction time 
8B; = y-intercept and coefficients to be determined 
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Figure 57: Mean Location Error And Standard Deviation (1987-1990) 
[Ref. 48] 
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could be found where the early TIP data could be used to approximate the final decay 

prediction time. [Ref. 48:pp. 22-23] The author states the following about equation 

(199): 
...multiple linear regression was then used to first determine if a model in the form 
shown...could be found to predict the final decay time. The results were an R- 
squared value of 1.0000 and a p-value of .0001. This means that at a significance 
level of .05 there exists a perfect linear relation between some of the independent 
variable and the dependent variable where at least two of the B terms are not zero. 
The variance inflation values were all extremely large, indicating the independent 
variables were all highly correlated and that a great deal of redundancy exists in 
the data. [Ref. 48:pp. 37-38] 

Based on the above results, six separate linear models were developed in an effort to 


approximate the final TIP with greater accuracy than the current process. Specific 
characteristics used in the modeling include the following : [Ref.48:p. 38] 
1. The first model uses only 7-day TIP data to calculate the expected value of the 
final decay prediction time, E(t). 


2. One additional decay prediction data point is incorporated in each subsequent 
model. 


3. All 180 TIP objects were used in the six models to calculate E(t). 


The six models are defined as 


Et) = -0.116442 + 0.999064(,) (200) 
Et) = -.0155478 + 0.49007(t,) + 0.951197(,) (201) 
E,t) = -.0082444 + 0.41692(t,) + 0.050001(t,) + 0.0908343(t,) (202) 


Pe )\| 


E,t\ = -0.053030 + 0.007471(¢,) - 0.014395(¢,) + 0.196518 (t,) 
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+ 0.810492(,) = 

Et) = -0.022322 - 0.004222(t,) = 0.005502(¢,) * 0.046458(r,) (204) 
-0.049000(,) + 1.012311() 

Egt) = -0.008370 + 0.000553 (¢,) - 0.002944(r,) + 0.01454 1(¢,) (205) 


-0.009096(t,) - 0.183413(t,) + 1.180378(,) 


where 
t, = 7-day prediction time 
t, = 4-day prediction time 
t; = 2-day prediction time 
t, | = 1-day prediction time 
ts = 12-hour prediction time 
ts = 6-hour prediction time 


Figure 58 shows the mean approximate error for the six regression models 
[Ref. 48:p. 40]. By comparing Figure 58 with the mean time error in Figure 56, the 
regression models show a better mean approximation error than the TIP runs. The 


conclusions and recommendations from the thesis are: [Ref. 48:p.42] 


1. The decay predictions were much better in general than the reported + 20 %. 
2. The use of linear models in conjunction with the data generated by the TIP 


processing would allow the SSC to better predict final decay time by elimination 
of positive bias in the data. 
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Figure 58: Regression Model Mean Approximate Error (1987-1990) 
[Ref. 48] 
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3. The SSC should conduct a study of the current Special Perturbations model to 
attempt to better account for the level of solar activity and its effect on the 
atmosphere. 

B. MONTE CARLO ANALYSIS OF SKYLAB’S IMPACT AREA 

Reference [45] investigated the debris impact point dispersion area of Skylab. A 
combination of Monte Carlo statistical analysis and parametric methods were used to 
determine the three-sigma limits of the debris footprint. The investigation was conducted 
by Martin Marietta Aerospace during the design and development of the Teleoperator 
Retrieval System (TRS) for the Skylab reboost/deboost mission. The dispersion analysis 
was conducted to support the deboosting of Skylab to a safe oceanic impact area clear 
of islands and routine shipping lanes. [Ref. 45:pp. 1-2] 

The Monte Carlo statistical analysis is an efficient and realistic approach that can 
be used to calculate the debris impact area because of the large number of input variables 
and nonlinearities of the problem. Impact dispersion area boundaries depend on the 


following factors: [Ref. 45:p. 2] 


1. Entry dispersions 

2. Relative flight path angle 
3. Debris ballistic coefficient 
4. Breakup altitude 


5. Environmental conditions (wind direction/magnitude and atmospheric density) 
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The entry dispersions were represented by a 6 X 6 covariance matnx of velocity 
and position errors at the entry altitude. For the deboost mission, the error sources were 


the uncertainties associated with: [Ref. 45:p. 2] 


1. Accelerometers 

2. Gyros 

3. Computer errors 

4. Initial alignment errors 


5. Vehicle performance dispersions 


The covariance matrix of state variables was determined by conducting a trajectory error 
analysis using a 1-sigma deviation of the errors and the nominal flight trajectory. In the 


application of the Monte Carlo analysis the major considerations were: [Ref. 45:p. 2] 


1. The simulation of perturbed trajectories. 
2. The modeling of the satellite breakup. 


3. The determination of impact points for each trajectory. 


At the entry altitude, the Monte Carlo method required generation of random state 
vectors from the 6 X 6 covariance matrix. Because the velocity and position errors are 
correlated, the random error vector in an uncorrelated space was derived from a 


transformation to a principle axes system. This vector is given by 


215 


e: = [nve,. Rees 5 Tel. (206) 


where 
nN; | =a set of random numbers drawn from a normal distribution of mean 0 and 
variance 1 
e; = the square roots of eigenvalues of the covariance matrix 


The original coordinate system random state vectors were derived from the equation 














R 
| se (207) 
pitty | tlellel 
Pp n 
where 
d = 6 X 6 eigenvector matrix that transforms the principle axes system into 
the original coordinate frame 
R, = perturbed radius vector 
V, = perturbed velocity vector 
R, = nominal radius vector 
V, = nominal velocity vector 


In order to calculate the dispersion area, the Monte Carlo analysis required a large 
number of simulated trajectories and random perturbed state vectors at the entry point 
for each perturbed initial condition. [Ref. 45:pp. 2-3] 

Additionally, the Monte Carlo scheme required a fast and efficient computational 
technique to calculate flight trajectories from several initial conditions. Analytical 
solutions, such as Loh’s second-order theory (presented in Chapter III) were not suitable 


for footprint dispersion analysis because these solutions are valid for only certain portions 
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of the trajectory and entry conditions. The trajectory simulation used the equations of 


motion of a nonlifting vehicle in the form 





7 pViV, + wy)g 
V, = 8, - ———— (208) 
a 2BP 
, pViV, + wx)g 
lao — (209) 
TS 2BP 
pViV sg 
V=¢9 - z (210) 
a 2BP 
where 
r= y(V, + wy)? + (Ysa wx) + (V,) (relative velocity) (211) 
2.02y82 = gravitational acceleration components 


V,,V,,V¥, = inertial velocity vector components 


Ky ,Z = radius vector components 
w = Earth’s rotation rate 
BP = ballistic parameter or coefficient (W/C,)A) 


Equations (208) through (210) were integrated several times from entry to the impact 
point for the Monte Carlo analysis. The integration step size used was as large as 
possible due to a stability consideration. The system approached a dynamic equilibrium 
condition which resulted in a numerical instability. This was a result of the opposing 
gravitational acceleration and aerodynamic deceleration during the vertical portion of the 


flight. A numerical search was conducted to determine the step size as a function of 


ZA) 


satellite debris ballistic coefficient. [Ref. 45:pp. 3-4] Figure 59 shows the stable and 
unstable regions used for the numerical integration [Ref. 45:p. 4]. 

Since the Monte Carlo analysis required the simulation of a large number of 
trajectories, an accurate computationally fast atmospheric density model was needed. 
The authors chose the 1962 U.S. Standard Atmosphere that provided density values from 
0 to 400,000 ft. [Ref. 45:p. 5] 

Although both static and dynamic global atmospheric models accounting for 
diurnal, seasonal, and latitudinal variations are available, it is found that such 
models are not suitable for the Monte Carlo dispersion analysis of satellite 
footprints. A sensitivity study of density variation indicates only a second order 
effect on footprint dispersion. Furthermore, the results of Purcell and Barbery 
show the downrange impact errors resulting from atmospheric variations are slight. 
[Ref. 45:p. 5] 

The simulation assumed one breakup altitude where all the pieces had the same 
initial position and velocity at breakup. By using a bounded parametric approach, the 
smallest and largest debris ballistic coefficients were used to determine the largest 
uprange and downrange footprint dispersions from the nominal impact point. Since the 
ballistic coefficient will vary as it passes through the atmosphere due to variations of C, 
and Mach number, the program used a generalized Cp) vs Mach number curve for 
tumbling pieces to update the ballistic coefficient at various altitudes. This curve is 
based on the results of a range safety study of Titan launch vehicle debris. [Ref. 45:p. 
5] 

For the footprint dispersion study, it is found that the transonic flow region, where 
Cy variation with Mach number is significant, occurs during the vertical descent 
of the debris with minimum effect on footprint dispersion. However, the variation 


of ballistic parameters with altitude affects the time of arrival of the debris on the 
ground. [Ref. 45:p. 5] 
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Figure 59: Numerical Stability Region 
[Ref. 45] 
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The latitude and longitude (lat/long) impact points of the smallest and largest 
Skylab debris were provided by the trajectory simulation for the various initial 
conditions at the entry altitude. The Lat/Long impact points were statistically treated 
which determined a three-sigma impact area boundary. The point estimate and 
confidence level estimate of the impact points for any given confidence level, e, and 
probability, a, are the fundamental statistical quantities of interest. 

These statistical computations relate the calculated longitude and latitude bounds 
from a finite sample to their true values corresponding to an infinite sample 


required for the Monte Carlo technique. [Ref. 45:p. 6] 


By using a normal distribution approximation given by 


2 
aK, (212) 
Cea) 


where K, 1s defined from the equation 





t 
ere 2 "5 (213) 
5 (1 - (2x) he e * dt 
equation (212) provides the estimated number of samples required to estimate a good 


point. A probabilistic statement given by 


P(O < 0.) = a (214) 
P(os,) = @ (215) 
where 

@ = impact point latitude 


6 impact point longitude 
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was obtained from the frequency histogram and cumulative probability distribution of 
Lat/Long impact points generated from the Monte Carlo analysis. For any specified 
value of a, the object of the statistical analysis was the estimation of ¢, and 6,. [Ref. 


45:p. 6] 


1. Simulation Results 
During the simulation to determine the debris footprint dispersion, Skylab was 

deboosted from a 170 nm circular orbit. The nominal entry and corresponding impact 
point Lat/Long were respectively (30°S,30°W) and (49.86°S,34.92°E). The nominal 
entry point state vector and the lower half of a symmetrical ECI frame covariance matrix 
are shown in Table XVI [Ref. 45:p. 7]. By using a Gaussian distribution for the Monte 
Carlo analysis, the initial state vectors were randomly selected from the covariance 
matrix. Specifically, 500 random entry states were created by using 3000 random 
numbers drawn from the distribution. 

To ensure that these numbers truly represent a normal distribution, their mean and 

variances were determined and adjusted to be 0 and 1 within stringent tolerances. 

The number of samples were found to be sufficient because further increase in 

sample size did not significantly effect the output variable distribution and 

probability limits. [Ref. 45:p. 8] 
Figure 60 shows a scattergram of the state vectors which define the entry flight path 
angle [Ref. 45:p. 7]. The trajectories were then simulated from the entry points down 
to the breakup altitude for each randomly selected state vector. Figures 61 and 62 show 


the results of a nominal run used to determine the flight characteristics of Skylab [Ref. 


45:pp. 8-9]. From Figure 62, the following parameters vs altitude were presented: 
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Table XVI: ENTRY STATE VECTOR / ERROR COVARIANCE MATRIX 
(Ref. 45] 


1 


Nominal Entry Point 


X= 17,648,032 ft V = -3238.58 ft/sec 
y= -5,174,810 ft ‘, = 21011.96 ft/sec 
z = -10,760,572 ft io -14693.59 ft/sec 
ECI Error Covariance Matrix 
X y Z My Vi Me 
x  5.76546E8 


-3.95374E9 2. /2016E10 
2.79469E9 -1.92075E10 1.35675E10 
4,86448E6 -3.34635E7 2.36298E7 4.167964 


—<=-— NO OMX 


V -],35955E6 9.34064E6 -6.59878E6 -1.14913E4  3.20968E3 
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Figure 60: Relative Entry Flight Path Angle Scattergram 
[Ref. 45] 
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Figure 61: Skylab Nominal Flight Characteristics 
[Ref. 45] 
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Figure 62: Downrange Impact Point Variations 
[Ref. 45] 
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In Figure 62, the 275,000 ft (84.8 km) breakup altitude of the simulation is shown along 
with the downrange variations for several debris sizes (denoted by ballistic parameter or 
coefficient). The lat/long histogram and cumulative probability distribution diagrams 
respectively shown in Figures 63 and 64 [Ref. 45:pp. 11-12] indicate the two distinct 
regions corresponding to the largest and smallest ballistic coefficient debris. 
Additionally, these diagrams determine the lat/long bounds where all debris is likely to 
fall. [Ref. 45:pp. 7-12] 

The downrange and crossrange debris impact dispersion distances were 
calculated using the lat/long and azimuth angles at the nominal impact point by the 


matrix equation [Ref. 45:p. 12] 
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Figure 63: Latitude Histogram / Cumulative Probability Distribution 
[Ref. 45) 
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Figure 64: Longitude Histogram / Cumulative Probability Distribution 
[Ref. 45] 
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| AE -CAS$,C8,-SAS8,  SACO,-CASH,S0,  CACh,| | Xp-Xy 
) 
| 


An | =| CASO, -SAS,CO,  -SASHSO0,-CACO,  SACH, ie be saa 
Ag Coh,CO,, Coh,S0, So, Zool, 
where 

En, ¢ = orthogonal coordinate system established at the nominal point 

E = axis tangent to the ground trace at the nominal point 

7 = axis along the radius vector 

c = axis completes a nght-handed triad 

on = nominal impact point latitude 

6 = nominal impact point longitude 

N = azimuth angle ( sin’[cos i / cos ¢y] ) 

Ay = RyConCOn 

Yn = RyCoySOy 

Zn = RySon 

Xp = RpCdpCé> 

Yp = RpC¢,S8, 

Zp = RpSdp 

Ss = sin 

c = Cos 


apis 


N nominal 


P = perturbed 
Figure 65 shows the downrange and crossrange dispersion distances from the nominal 
impact point calculated from equation (216) [Ref. 45:p.14]. As shown in the figure, the 
crossrange dispersions were relatively small as compared to the downrange dispersions. 
In order to determine the effect of the breakup altitude on the dispersion area, 
Monte Carlo simulations were conducted at various altitudes from 200,000 to 350,000 
ft (61.7-107.9 km). In this altitude range, the reentry object experiences large thermal 
and structural loads. Figure 66 shows three-sigma downrange and uprange dispersion 
for the largest and smallest debris [Ref. 45:p. 15]. 
This figure is useful for estimating the total footprint dispersion resulting from 
breakup of smaller pieces at higher altitudes and heavier pieces at lower altitudes, 
thereby accounting for the uncertainty in breakup altitude. [Ref. 45:p.14] 
In their conclusions, the authors note significant footprint dispersion 
variations for breakup in the 200,000 to 350,000 ft altitude region. Furthermore, they 
state that the comprehensive Monte Carlo analysis is very useful and appropriate in 


determining impact dispersion areas of a spacecraft or discarded portions of a launch 


vehicle. [Ref. 45:p. 15] 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The primary goal of this thesis was to identify the "state-of-the-art" of orbit-decay- 
induced uncontrolled reentry and impact prediction methods, with an emphasis on the 
physics of the final few revolutions to impact. This was accomplished through a 
comprehensive literature survey from the 1950’s to the present of unclassified military 
and civil databases. The survey indicated that there is some significant foreign work 
being done and much of it was not available to the authors, in English translation, and 
thus was not included in this survey. Also, the literature survey reflects the changing 
scientific terminology over the course of several decades and it is especially noticeable 
in the different forms that the common variables take in the numerous equations 
presented. The authors did not make any attempt to use the standard AIAA 
astrodynamics nomenclature or to standardize the equations in any other way. 

The principal conclusion of this thesis is that the current uncontrolled reentry and 
impact prediction methodology, used in the U.S. and abroad, is based on analysis which 
is 30 or more years old. This conclusion is based on the fact that the U.S. method takes 
its roots in the works of Brouwer (1959) and Allen and Eggers (1957), and that the U.S. 
method is the accepted international standard, as shown by the literature survey. 

While conducting the literature search dating back to the 1950’s, the authors 


noticed a definite trend, through the years, in the focus of published material related to 
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reentry. During the timeframe from the late 1950’s until the mid 1960’s, the emphasis 
in the literature was primarily devoted to understanding and describing the physics of 
reentry, specifically as it pertains to controlled ballistic reentry of Mercury and missile 
type vehicles. With the exception of Sputnik IV, very little of the literature surveyed 
during this timeframe was strictly devoted to uncontrolled reentry. However, during the 
derivations of their analytical reentry theories, some early pioneers such as Chapman and 
Loh investigated the shallow orbit-decay-induced reentry. With the physics of reentry 
fairly well understood, starting in the late 1960’s and continuing into the 1980’s, the 
emphasis in the literature had shifted to controlled, gliding reentry in order to support 
the launch of the Space Shuttle. In 1965, Eggers and Cohen stated: 
Significant advances have been accomplished in the science and technology 
appropriate to atmosphere entry of spacecraft during the eight years since the 
launching of Sputnik I. This progress is illustrated by the successful entry from 
Earth orbit of the manned Vostok, Mercury, Voskhod, and Gemini spacecraft, for 
which the problems in orbital entry, such as high convective heating rate and load 
and communication blackout, were successfully over come...Although the first 
Apollo vehicle entry has yet to be demonstrated... Much engineering work for this 
vehicle remains to be done; however, the fundamental research activity associated 
with Apollo is being reduced in favor of that associated with missions and vehicles 
of the more distant future. [Ref. 75:pp. 339-340] 
Uncontrolled reentry briefly came back into focus in the late 1970’s and early 1980s 
with the reentry of Cosmos-954 and Skylab and continued throughout the 1980’s in 
varying degrees where the reentries of Cosmos-1402, 1601 and Salyut-7/Cosmos-1686 
continued to spark a flurry of literature from the European Space Agency. 


The literature survey of recent publications shows a strong reliance on work done 


in the 1950’s and 1960’s as a basis for "extensions" or "modifications" of pre-existing 
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methods used for prediction of reentry and impact. Although reliance on pre-existing 
methods is not in and of itself flawed, the modern works often subscribe to the same 
assumptions which the original authors were forced to make because of hardware or data 
limitations of their time. This serves to perpetuate the inherent limitations of the 
method’s ranges of applicability and validity. For example, it has been shown [Refs. 
14,16-17] that some very important work done in the 1950’s, such as Allen and Eggers 
[Ref. 65], which is still routinely referenced in recent publications, is inherently flawed 
due to assumptions made under conditions of little or no data. This particular reference 
is especially important since it is the basis for most of the modern reentry heating work. 
It becomes even more significant since the literature survey has shown the coupled 
dependence of reentry breakup to reentry heating and dynamic load effects. Reference 
[42] describes the limitations imposed by these assumptions through mathematical proofs. 
This is further substantiated by reentry breakup observations [Ref. 86] which indicate that 
the classical convective heat transfer equations, when applied to breakup analysis, 
consistently underestimate reentry survivability. 

Of the various "extensions" to the current reentry theory, of which the NORAD 
method is recognized as the international standard, there does not appear to be any one 
method which is singularly superior to the others. However, it is the opinion of the 
authors that the ESA FOCUS program merits special attention for further research. This 
conclusion is based on the fact that this program contains a very sophisticated method for 
dealing with the drag coefficient, C,), and could easily incorporate other higher-order 


effects occurring in reentry, such as attitude generated lift, rotation of the orbital plane 
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(due to atmosphere rotation at attitudes < 60 km), linear variations in the ballistic 
coefficient in the final few days of orbit and others as mentioned previously in this thesis. 

Numerous works cited in the literature survey also describe computing power and 
computer time constraints as critical parameters of the problem. These hardware 
limitations then forced one or more simplifying assumptions to be made in order to deal 
with the technology limitations of the day. Subsequently, much of the literature survey 
shows work which has limited applicability for the problem of concern. 

Additionally, there is a general lack of consensus in the literature as it pertains to 
orbit-decay-induced uncontrolled reentry. The very definition of "reentry" remains vague 
and often is defined differently by investigators attempting to do similar research. Many 
investigators have examined various aspects of the problem; however, when surveying 
the literature, a common approach and standardized starting point for solving this very 
dynamic problem is not apparent. A good example to illustrate this point is the effect 
of angle of attack, a, on the ability to accurately predict uncontrolled reentry and impact. 
The ability to properly model the uncertainties in the reentry body configuration such as 
changes in area, mass, and attitude is a very difficult problem. As a result of these 
changes, three-dimensional angle of attack becomes an exceedingly difficult parameter 
to characterize especially when the reentry body is undergoing rapid configuration 
changes due to ablation and structural deformation. However, in order to properly model 
the aerodynamic coefficients for lift and drag, angle of attack must be considered. 
Specifically, lift as a function of angle of attack not only affects the trajectory (location 


of impact), but it also affects the altitude of breakup (survivability and dispersion area) 
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since it reduces the thermal and structural loads on the vehicle. However, in most cases 
where angle of attack is investigated, this coupled effect is not mentioned nor pursued. 
Even more likely in the literature, the reentry body is modeled as a point mass where 
angle of attack is not considered or lift is assumed to be zero or negligible. 

Another observation from the survey is that no comprehensive sensitivity and error 
analysis has been conducted in order to determine quantitatively the effects of critical 
parameters/variables on impact prediction. For example, numerous sensitivity and error 
analysis studies have been conducted on atmospheric density models uncoupled from 
other critical parameters. However, a coupled analysis of all of the critical parameters 
could determine which variables contribute most significantly to the overall accuracy of 
the various impact prediction methods. 

In addition to the author’s observations and conclusions described in the previous 
paragraphs, a summary of specific conclusions from the literature survey are: 

1. There is no clear definition of when reentry occurs. Typically 120 km or an 
orbital period of 87.5 minutes or less is used in the various models to define the 
Start of reentry. However, observations and post flight analysis indicate a range 
of altitudes where "reentry" actually occurs. 

2. There is a lack of observation data in the reentry regime due to a lack of global 
sensor coverage. This lack of data significantly adds to the uncertainty of the 


problem. 


3. The current deterministic dynamics model appears to be inadequate for 
processing the true physics of reentry. 


4. The lack of observational data coupled with the inadequate knowledge of the 


true physics of uncontrolled reentry significantly increases the uncertainty of the 
problem. 
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5. Reentry is most likely to occur within several tenths of a degree of the minimum 
altitude (height above the ground) point in the orbit, which i is not necessarily the 
point of perigee. [Ref. 99] 

6. The uncontrolled reentry is three times more likely to occur in a latitude band 
equaling the inclination than in an equatorial band (0 to 10 degrees), with the 
exception of polar orbits where an equatorial reentry is most likely to occur. 
[Ref. 99] 

7. The impact location (downrange and crosstrack) is affected by: [Refs. 45,101] 
(a) the rotation of the atmosphere, starting at an altitude of 30-60 km 
(b) the debris ballistic coefficient 
(c) the entry dispersions (velocity and position errors at the entry altitude) 
(d) the relative flight path angle 
(e) the breakup altitude 


(f) the environmental conditions such as wind direction/magnitude and 
atmospheric density 


8. There appears to be linear variations in the ballistic coefficient in the final few 
days prior to reentry. [Ref. 92] 


B. RECOMMENDATIONS 
Based on the above conclusions, the following recommendations are proposed: 


1. As mentioned in Chapter I, there was no attempt made to standardize the 
nomenclature and variables of the equations presented in this thesis, in 
accordance with AIAA astrodynamics standards or any others. Based on the 
literature survey, there clearly exists a need for such standards, especially when 
dealing with work spanning the course of several decades. 
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This literature survey was conducted through unclassified sources, written in 
English, only. A survey of classified and foreign language sources should be 
conducted. 


The primary focus of this thesis was a thorough survey of the physics of 
uncontrolled reentry. Throughout the survey, the authors noted several 
operational issues which should be addressed in follow-on research. These issues 
include human factors (orbit analyst experience), sensor bias and coverage, and 
improved international information exchanges. 


Several reentry/impact prediction methods were presented (LIFETIME, FOCUS, 
and NORAD) in this thesis. A side by side comparison of these methods using 
historical data should be conducted in order to compare and contrast their 
overall accuracies and efficiencies. 


A sensitivity and error analysis of critical parameters including atmospheric 
density, aerodynamic coefficients, initial conditions (vehicle area, mass, attitude, 
position and velocity), thermal and dynamic loads, and breakup phenomena 
should be conducted in order to improve the current dynamics model and to 
focus future research efforts. 


In view of the lack of global sensor coverage, especially in the southern 
hemisphere, high priority reentry events may be better predicted through a joint 
cooperative effort using Navy and Air Force assets as mobile supplemental 
observation and tracking stations. 


Since uncontrolled reentry is characterized by numerous rapidly changing 
critical parameters, which are ill-defined, stochastic and statistical methods 
should be applied to current reentry models to better analyze the sensitivity of 
the various uncertainties associated with this problem. This could serve as an 
“operational tool" to help improve the prediction accuracies of the current 
models while the physics of the reentry process is being more thoroughly 
investigated. 


The authors strongly recommend future cooperative research between the Naval 


Postgraduate School and the AFSPACECOM/USSPACECOM in order to solve 
this largely unknown and highly dynamic process. 
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